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1.  Introduction 


This  report  describes  recent  advances  in  the  analysis  of  radar  target  signatures 
for  the  purpose  of  target  identification.  The  goal  of  the  research  has  been  to  de¬ 
velop  techniques  for  automatic  classification  of  radar  targets  by  digital  processing  of 
measured  scattering  signals. 

Many  approaches  have  been  considered  for  target  identification  systems.  Of  these 
approaches,  we  have  concentrated  on  processing  the  scattering  data  to  extract  time 
(range)  domain  scattering  features.  In  particular,  we  are  interested  in  the  ranges  of 
the  dominant  scattering  centers,  and  the  amplitude  (and  polarization)  properties  of 
each  scattering  center.  These  scattering  center  ranges,  along  with  their  associated 
properties,  would  form  a  feature  vector  for  use  in  a  target  classification  system. 

There  are  two  main  approaches  for  extracting  the  ranges  of  scattering  centers 
from  measured  frequency-domain  data.  Traditionally,  target  scattering  centers  are 
found  by  using  Fourier  transform  techniques  to  convert  the  frequency-domain  data 
into  a  time-domain,  or  range-domain,  profile  [l].  This  profile  is  then  searched  for 
sharp  peaks,  and  the  ranges  of  these  peaks  form  the  scattering  center  range  esti¬ 
mates.  The  amplitudes  of  the  peaks  give  the  magnitudes  of  the  scattered  field  for 
each  scattering  center  when  single  polarization  data  is  available.  Little  work  has  ap¬ 
peared  on  Fourier  processing  of  multiple  polarization  data;  one  method  for  extracting 
polarimetric  properties  of  scattering  centers  was  recently  reported  in  [2]. 

This  report  considers  a  different  approach  for  scattering  center  extraction:  the 
parametric  modeling  approach.  In  parametric  modeling,  the  target  scattering  char¬ 
acteristics  are  described  by  a  model,  and  the  parameters  of  that  model  are  then 
estimated  from  measured  data.  The  model  is  chosen  so  that  the  model  parameters 
relate  to  target  scattering  centers  in  a  natural  way.  In  addition,  the  model  we  have 
chosen  extends  naturally  to  the  case  in  which  data  for  multiple  polarizations  are 
available. 
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An  outline  of  the  report  is  as  follows.  Section  2  describes  the  parametric  modeling 
approach  in  detail,  and  discusses  the  advantages  of  such  an  approach.  Section  3 
outlines  the  technical  accomplishments  of  this  research  effort.  Section  4  concludes 
the  report. 
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2.  Technical  Approach 


Assume  that  a  given  radar  acquires  stepped  frequency  measurements  of  the  (co¬ 
herent)  scattering  from  a  target.  Denote  these  measurements  as 

sxy{fk),  k  =  0,1,. . .  ,N  -  l.  (2.1) 

Here,  the  subscripts  x  and  y  denote  the  received  and  transmitted  polarizations  of 
the  radar,  respectively;  //,  denotes  the  frequency  at  the  fcth  frequency  step.  For  the 
purpose  of  clarity,  we  assume  that  the  radar  transmits  left-circularly  polarized  signals 
and  receives  both  horizontally  and  vertically  polarized  fields;  thus,  both  Sh.i(fk)  and 
svi{fk)  are  measured.  However,  other  polarizations  could  be  used. 

The  parametric  modeling  method  involves  modeling  the  target  as  a  set  of  scat¬ 
tering  centers.  In  particular,  we  use  the  Prony,  or  damped  exponential  model  [3,  4]. 
In  the  frequency  domain,  the  model  takes  the  form  of  sum  of  damped  exponential 
terms: 


i>u(Jk) 

II 

Ah.li 

Pi*, 

0  <  k  <  N  -  1 . 

(2.2) 

!=1 

Ayh 

Here,  the  A ,  and  p,  parameters  represent  the  amplitude,  range,  and  damping  factor  of 
each  of  the  M  scattering  centers.  Each  scattering  center  is  characterized  by  a  pole  p,-, 
whose  angle  Zp,  gives  the  range  of  the  ith  scattering  center,  and  whose  amplitude  |p,| 
relates  to  the  range  dispersion  of  the  scattering  center.  For  an  ideal  point  scatterer, 
|p, |  =  1,  but  for  realistic  targets  the  scattering  will  be  attenuated  slightly  as  frequency 
either  increases  or  decreases,  thus  |p,  |  will  vary  a  little  bit  around  one.  The  complex 
amplitudes  A^it  and  Avu  characterize  the  amplitude  and  polarization  properties  of 
the  ith  scattering  center;  the  amplitude,  ellipticity,  and  tilt  of  the  scattering  center 
can  be  found  by  a  simple  transformation  of  Aui  and  Avh  [3].  This  model  has  the 
advantage  that  full  polarization  information  is  incorporated  into  a  single  model;  on 
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the  other  hand,  if  not  all  of  the  polarization  data  is  available,  the  reduced  data  set  can 
be  used  to  obtain  an  exponential  model  with  fewer  attributes  per  scattering  center. 

Parametric  modeling  of  target  signatures  offers  several  advantages  over  traditional 
Fourier-based  target  signature  analysis.  These  include: 

•  High  Resolution:  Parametric  methods  offer  higher  resolution  than  can  be 
obtained  by  Fourier  processing;  as  a  result,  lower  bandwidth  radars  may  be 
used  to  achieve  a  desired  range  resolution. 

•  Automatic  Feature  Extraction:  The  parametric  methods  we  propose  di¬ 
rectly  estimate  the  location  of  scattering  centers,  and  directly  estimate  a  set 
of  parameters  which  characterize  each  scattering  center.  This  is  in  contract  to 
Fourier-based  processing  techniques,  in  which  some  sort  of  feature  extraction 
procedure  must  be  applied  to  the  estimated  range  profile  of  the  target. 

•  Separate  Characterization  of  Each  Scattering  Center:  The  proposed 
parametric  models  characterize  each  scatterer  with  a  set  of  attributes  which 
include  polarization  properties,  frequency  response,  and  in  some  cases  a  de¬ 
scription  of  the  physical  nature  of  that  scatterer. 
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3.  Technical  Accomplishments 


Below  we  outline  the  technical  accomplishments  developed  under  this  research 
program. 

3.1  Statistical  Analysis  of  TLS-Prony  Techniques 

The  main  algorithm  we  used  for  estimating  the  parameters  of  the  exponential 
model  is  the  so-called  TLS-Prony  estimator.  This  method  involves  setting  up  a  set  of 
backward  linear  equations,  and  solving  for  prediction  coefficients  using  a  total-least- 
squares  (TLS)  techniques.  The  TLS  solution  involves  performing  a  singular  value 
decomposition  on  a  data  matrix.  The  zeros  of  the  prediction  polynomial  associated 
with  the  backward  linear  prediction  coefficients  form  the  exponential  pole  estimates. 
Once  the  exponential  terms  are  found,  the  corresponding  amplitude  coefficients  are 
then  determined  by  a  least-squares  technique.  This  method  is  attractive  because  it  is 
computationally  simple  (involving  linear  algebraic  solutions  of  equations  and  a  single 
polynomial  rootfinding  step),  but  gives  estimates  whose  statistical  properties  are  close 
to  optimal. 

Our  first  main  contribution  is  a  complete  statistical  analysis  of  the  poles  and 
amplitude  coefficients  obtained  lrom  the  TLS-Prony  technique  when  the  measurement 
data  are  corrupted  by  Gaussian  white  noise.  The  analysis  gives  the  statistics  of 
the  model  coefficients  for  high  SNR.  Because  the  pole  and  amplitude  parameters 
relate  directly  to  physical  properties  of  the  scattering  centers,  this  analysis  gives  a 
theoretical  prediction  of  scattering  center  extraction  performance  in  the  presence  of 
noise  (see  also  [5]).  In  particular,  resolution  of  scattering  centers,  and  accuracy  of 
scattering  amplitudes  and  polarization  properties  can  be  found  as  a  function  of  SNR 
and/or  polarization  diversity.  The  statistical  theory  is  compared  to  Monte-Carlo 
simulation  results,  and  lound  to  agree  quite  well  even  for  moderate  SNRs  (-5  dB  to 
5  dB  depending  on  the  signal  scenario).  These  results  have  been  published  as  a  M.S. 
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Thesis  by  Mr.  Ching-Hui  Ying  [6],  and  also  as  a  conference  paper  presented  at  the 
1991  Asilomar  Conference  [7].  These  results  have  also  been  accepted  for  publication 
in  archival  Journal  Automatica  in  a  special  issue  on  Statistical  Signal  Processing  and 
Control  [8].  The  Automatica  paper  is  included  in  Appendix  A. 

3.2  Computationally  Efficient  TLS-Prony  Estimation  Algorithms 

A  second  advance  is  the  development  and  analysis  of  TLS-Prony  estimation  algo¬ 
rithms  which  have  lower  computational  cost  then  the  standard  TLS-Prony  method, 
but  which  have  almost  the  same  accuracy  performance.  Specifically,  the  method 
involves  decimation  of  the  frequency-domain  scattering  measurements.  Decimation 
reduces  the  length  of  the  data,  and  gives  a  large  corresponding  reduction  in  compu¬ 
tational  cost. 

Specifically,  the  data  is  processed  through  several  different  filtering  and  decima¬ 
tion  blocks.  Each  of  these  blocks  is  designed  to  estimate  the  scattering  centers  of 
the  target  for  a  particular  segment  of  the  range  profile.  The  individual  segments 
can  be  processed  in  parallel,  resulting  in  a  significant  savings  of  computation  (at  an 
increase  of  processing  hardware).  Alternately,  the  processing  can  be  done  sequen¬ 
tially.  In  either  case,  the  total  number  of  computations  is  substantially  less  than  the 
computations  involved  for  a  single  estimate  of  all  scattering  centers. 

One  concern  with  such  an  approach  is  that  the  statistical  accuracy  of  the  resulting 
scattering  center  estimates  is  lower  than  the  accuracy  obtained  in  a  single  estimate 
approach.  To  address  this  issue,  we  derived  analytical  expressions  which  predict  the 
accuracy  of  the  scattering  center  range  and  amplitude  (and,  in  the  case  of  multiple 
polarization  data,  the  polarization  characteristics  of  the  scattering  center  as  well) 
for  the  decimation-based  TLS-Prony  algorithm.  This  analysis  again  assumes  high 
SNR,  but  is  shown  through  comparisons  with  Monte-Carlo  experiments  to  hold  for 
moderate  and  low  SNRs  as  well.  The  analysis  and  simulation  demonstrate  that 
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the  decimation-based  algorithm  can  provide  scattering  center  estimates  with  almost 
the  same  accuracy  as  in  the  single  estimation  algorithm  (the  loss  is  about  2-4  dB), 
but  at  huge  computational  savings  (reductions  by  a  factor  of  20,  using  sequential 
computation  on  a  single  CPU,  were  for  typical  examples;  if  parallel  computation 
hardware  is  used,  the  savings  would  be  greater  by  an  additional  factor  of  about  6- 
10). 

The  decimation-based  algorithms,  and  the  statistical  analysis  of  these  algorithms, 
have  been  published  in  a  conference  paper  which  appeared  in  the  1992  SPIE  Aerospace 
Sensing  Symposium  [9],  and  has  also  been  submitted  to  the  IEEE  Transactions  on 
Signal  Processing  [10].  The  IEEE  Transactions  paper  is  included  in  Appendix  B. 

3.3  Detection  and  Resolution  of  Scattering  Centers 

This  research  is  concerned  with  detection  and  resolution  of  scattering  centers 
which  are  estimated  from  noisy  stepped-frequency  measurements  of  that  target.  The 
research  is  aimed  at  answering  the  following  questions: 

1.  As  a  function  of  radar  bandwidth,  SNR,  and  polarization  diversity,  how  closely 
spaced  can  scattering  centers  be  before  they  can  no  longer  be  resolved? 

2.  What  is  the  probability  of  detection  versus  probability  of  false  alarm  for  a 
scattering  center  (as  a  function  of  bandwidth,  SNR,  and  polarization  diversity)? 

The  question  of  resolution  is  closely  related  to  the  statistical  analysis  of  the  TLS- 
Prony  methods  as  described  previously.  Basically,  the  resolution  question  is  answered 
by  appropriately  applying  these  analysis  results  to  the  radar  scattering  application. 

The  detection  question  arises  from  the  mechanics  of  the  TLS-Prony  methods. 
The  TLS-Prony  methods  achieve  high  resolution  by  first  finding  a  backward  linear 
prediction  filter  of  high  order.  The  zeros  of  this  prediction  filter  are  found;  some  of 
these  zeros  correspond  to  scattering  centers  on  the  target,  and  some  are  extraneous, 
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“computation”  zeros.  One  must  then  decide  which  zeros  correspond  to  scattering 
centers;  this  decision  is  based  on  the  energy  of  the  “mode”  associated  with  that  zero. 
Because  the  data  are  noisy,  there  will  be  some  times  in  which  a  computation  zero  has 
energy  which  exceeds  the  decision  threshold,  giving  rise  to  a  false  alarm  scattering 
center.  Similarly,  there  will  be  times  in  which  a  true  scattering  center  has  energy 
which  lies  below  the  threshold,  and  the  scattering  center  will  not  be  detected.  By 
varying  the  threshold,  the  detection  versus  false  alarm  probability  can  be  determined 
for  various  radar  parameter  settings  (such  as  SNR,  polarization,  and  bandwidth). 

The  detection  question  cannot  be  answered  directly  in  terms  of  our  previous  analy¬ 
ses,  because  these  analyses  did  not  include  the  statistical  properties  of  the  extraneous 
poles.  Thus,  we  first  derived  the  statistical  properties  of  the  extraneous  poles;  this 
derivation  is  presented  in  [11].  From  this  analysis,  and  using  properties  of  the  en¬ 
ergy  of  a  scattering  center,  the  desired  detection  performance  measures  could  then  be 
calculated. 

A  conference  paper  which  describes  this  work  was  presented  at  the  International 
Radar  ’92  conference,  and  appears  in  the  conference  proceedings  [12].  The  This  paper 
is  included  in  Appendix  C. 

3.4  Scattering  Center  Extraction  from  Measured  Radar  Data 

Finally,  we  applied  our  parametric  modeling  results  to  some  measured  S-band 
Linear  FM  scattering  data  which  was  supplied  to  us  by  Rome  Laboratories.  Unfor¬ 
tunately,  only  a  limited  amount  of  data  was  available,  and  this  data  was  taken  at 
a  single  polarization,  so  we  were  unable  to  make  extensive  comparisons  between  the 
theoretical  statistical  analyses  and  the  results  obtained  from  measured  radar  data.  In 
addition,  we  had  no  calibration  standard  for  the  data,  so  it  was  not  possible  to  com¬ 
pare  the  estimated  scattering  centers  with  a  “ground  truth”.  However,  we  were  able 
to  conduct  limited  studies  on  scattering  center  accuracy,  and  on  bandwidth  extrapo- 
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lation  capability,  of  the  TLS-Prony  algorithms.  Based  on  these  studies,  we  found  that 
scattering  center  resolutions  between  2-4  times  the  resolution  capability  of  Fourier 
techniques  could  be  obtained. 

1  he  experimental  study  on  measured  radar  data  is  presented  in  detail  in  the 
technical  report  “Prony  Modeling  of  Linear  FM  Radar  Data,”  [11];  this  report  is 
included  in  Appendix  D. 
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4.  Conclusions 


This  report  presents  the  technical  advances  in  radar  scattering  center  extraction 
which  were  developed  under  this  contract.  We  have  made  significant  advances  in 
the  development  and  analysis  of  Prony-based  scattering  center  estimation.  We  have 
developed  theoretical  analyses  which  predict  the  performance  of  these  algorithms  for 
scattering  center  resolution,  accuracy  of  scattering  center  amplitude  and  polarization 
characteristics,  and  detection  versus  false  alarm  probabilities  of  scattering  centers. 
This  analysis  provides  these  performance  measures  as  functions  of  SNR,  data  length 
(which  is  directly  related  to  data  bandwidth),  and  data  diversity  (which  is  directly 
related  to  polarization  diversity  of  the  data).  The  analysis  assumes  high  SNR,  but  is 
shown  by  Monte-Carlo  simulations  to  be  applicable  for  low  to  moderate  SNRs  as  well. 
We  have  also  applied  these  algorithms  to  a  limited  set  of  measured  radar  data,  and 
shown  scattering  center  accuracies  of  two  to  four  times  that  obtainable  with  Fourier- 
based  techniques.  The  results  of  this  research  have  been  published  in  two  technical 
reports,  three  conference  papers,  and  two  archival  journal  papers. 

Future  work  should  focus  on  more  extensive  testing  of  these  methods  on  measured 
data,  especially  fully  polarimetric  data.  In  addition,  future  signal  processing  research 
on  constrained  TLS-Prony  or  other  algorithms  (in  which,  for  example,  the  poles  are 
constrained  to  he  on  or  near  the  unit  circle),  and  on  automatic  model  order  selection 
are  of  interest. 
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I.  Introduction 


The  problem  of  estimating  model  parameters  of  noisy  exponential  signals  is  an 
active  area  of  research.  This  problem  has  applications  in  a  number  of  areas,  including 
speech  processing,  deconvolution,  radar  and  sonar  signal  processing,  array  processing, 
and  spectrum  estimation.  A  number  of  authors  have  considered  various  aspects  of 
this  problem  [1]— [22] ,  and  a  large  number  of  algorithms  have  been  developed  and 
analyzed  [8.  11,  12],  [16]-[18],  [20,  21]. 

One  popular  class  of  algorithms  for  estimating  parameters  from  noisy  exponen¬ 
tial  sequences  are  the  subspace- based  approaches  [l]— [10],  [13]— [16],  [19]— [22] .  These 
include  the  MUSIC  algorithm  and  its  enhancements  [1,  16,  21],  subspace  rotation 
methods  such  as  ESPRIT  [5,  14,  15,  20,  21],  iterative  maximum  likelihood  meth¬ 
ods  [6,  7,  13],  minimum  norm  methods  [8,  19],  and  total  least  squares  (TLS)  meth¬ 
ods  [2,  9,  10,  22].  These  methods  have  proven  attractive  because  they  exhibit  good  sta¬ 
tistical  performance  at  a  modest  computational  cost.  This  has  been  well-established 
by  a  large  number  of  numerical  studies. 

More  recently,  there  has  been  interest  in  quantitatively  evaluating  these  methods. 
To  this  end,  a  number  of  researchers  have  analyzed  the  statistical  properties  of  such 
algorithms  [3.  8.  11,  12.  1C.  17.  18,  20,  21],  In  [3],  Henderson  presents  a  geometric 
study  ol  the  pole  estimation  problem,  and  analyzes  the  statistical  properties  of  the 
prediction  coefficients  when  the  Hankel  data  matrix  is  corrupted  by  an  i.i.d.  noise 
matrix.  Several  authors  have  presented  results  relating  to  pole  angle  (frequency) 
estimates  from  arrays  when  the  exponential  signals  are  undamped,  e.g.,  [8,  11,  16, 
17,  20,  21,  23].  A  related  perturbation  analysis  of  SVD-based  methods  is  presented 
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in  [24,  25,  26,  27]  and  applied  to  both  frequency  estimation  and  threshold  analysis 
for  exponential  modes.  In  [12]  Hua  and  Sarkar  present  an  angle-only  analysis  for 
the  least  squares  Prony  method  for  the  poles  of  undamped  exponentials.  Less  has 
appeared  which  considers  the  statistical  properties  of  the  parameters  for  damped 
exponentials.  Porat  and  Friedlander  [28]  consider  the  related  problem  of  ARMA 
system  identification  using  SVD-based  approaches.  Hua  and  Sarkar  [18]  present  an 
analysis  for  the  pole  estimates  of  damped  exponential  signals  using  their  matrix  pencil 
method,  but  have  not  presented  the  statistical  properties  of  the  amplitude  coefficients. 

This  paper  presents  an  extension  of  the  above  works  to  treat  a  general  exponential 
case.  We  introduce  a  complete  statistical  derivation  for  the  poles  and  amplitude 
coefficients  estimated  using  a  TLS-Prony  scheme  where  signals  consist  of  arbitrary 
damped  exponential  terms  in  noise.  We  provide  complete  statistics  for  the  individual 
pole  parameters  for  an  exponential  model  in  which  the  poles  may  lie  on,  inside, 
or  outside  the  unit  circle.  In  addition,  we  derive  the  statistical  properties  of  the 
amplitude  coefficients  associated  with  these  exponential  modes. 

The  Jesuits  of  this  paper  provide  a  sound  basis  for  performance  analyses  of  the 
TLS-Prony  estimation  method.  We  extend  previous  works  by  considering  the  general 
damped  case,  as  well  as  by  including  amplitude  coefficient  parameters  in  addition  to 
pole  parameters.  These  results  provide  the  tools  to  analyze  various  situations  and 
evaluate  the  potential  success  of  applying  the  TLS-Prony  estimation  algorithm. 

The  TLS-Prony  estimation  procedui’e  is  a  multi-snapshot  extension  of  the  algo¬ 
rithm  presented  in  [9,  10].  The  advantage  of  singular  value  decomposition  (SVD)  in 
noise  cleaning  of  the  Toeplitz  data  matrix  is  well-known.  The  multiple  snapshot  incor- 
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poration  is  a  straightforward  one  in  which  more  than  one  set  of  amplitude  coefficients 
corresponds  to  the  set  of  poles.  The  procedure  is  discussed  in  Section  II. 

The  statistical  derivation  for  this  procedure  is  based  on  a  first  order  perturbation 
analysis;  thus  the  analysis  assumes  high  SNR.  We  derive  the  complete  covariance  ma¬ 
trix  of  the  estimated  parameters  for  this  case.  The  parameters  include  the  magnitudes 
and  angles  of  the  poles,  and  the  magnitudes  and  angles  of  the  amplitude  coefficients. 

Using  these  expressions,  several  general  properties  of  the  parameter  covariance 
matrix  are  derived  for  the  high  SNR  case.  We  show  that  the  angle  and  magnitude 
parameters  are  uncorrelated  for  both  the  poles  and  the  amplitude  coefficients.  We 
also  show  that  if  the  relative  magnitude  of  the  pole  or  amplitude  coefficient  estimate 
is  considered  (i.e.  where  a  is  the  true  magnitude),  then  the  corresponding  angle 
and  relative  magnitude  variances  are  equal. 

This  paper  also  examines  pole  estimation  accuracy  as  functions  of  pole  magnitude, 
data  length,  and  pole  separation  using  the  variance  expressions.  We  compare  these 
variance  results  to  the  corresponding  Cramer-Rao  bounds  and  verify  the  theoretical 
results  using  Monte-Carlo  simulations.  The  effects  on  poles  inside  and  outside  the 
unit  circle  using  backward  or  forward  linear  prediction  in  the  TLS-Prony  estimation 
scheme  are  also  detailed. 

An  outline  of  this  paper  is  as  follows.  Section  II  presents  the  data  model.  Sec¬ 
tion  III  presents  the  statistics  of  the  model  parameters  and  their  properties.  Sec¬ 
tion  IV  presents  some  examples  using  the  statistical  expressions.  Finally,  Section  V 
concludes  the  paper. 
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II.  Estimation  Procedure  Review 


A.  Data  Model 

Assume  we  have  N  “snapshots”  of  data  vectors  y(t),  each  of  length  m: 

y(t) =  yo(t)  yi{t)  ■■■  ym-i(t)  t  =  l,2,...,N.  (1) 

Each  data  vector  is  modeled  as  a  noisy  exponential  sequence 

n 

y9(0  =  52xi{t)Pi  +  e,(0  q  =  o,  1.  (2) 

»=i 

There  are  n  distinct  exponential  modes  in  the  data;  the  n  poles  do  not  vary 

from  snapshot  to  snapshot,  but  the  amplitudes  x,(t)  may  vary.  Here,  it  is  assumed 
that  {eq(t)}  are  uncorrelated  zero  mean  complex  white  Gaussian  noise  sequences  with 
variance  a.  Equation  (2)  may  be  compactly  written  as 


y(t)  =  Ax(t)  +  e(t),  (3) 


r 

T 

1  T 

where  t (/.)  =  to(<)  <-,(/)  ■ 

■  em_,(0  ,  At)  = 

Xi{t)  x2{t)  ■■ 

••  arn(0 

and  A  is  the  m  x  u  Vandermonde  matrix  derived  from  n  signal  poles 
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B.  Parameter  Estimation 

The  multi-snapshot  backward  linear  prediction  equations  are  given  by: 


where 


and  where 


r  1 

y  ■  Y  ~0, 

J  b 


b  bi  62  •  •  •  bi 


y  :  V  = 


2/o  (H 

2/1  (1) 

2/2(1) 

2/i(l) 

2/1  ( 1 ) 

2/2(1) 

2/3(1) 

2/L  +  l(l) 

2/m-(t  +  l)(l) 

2/m-i(l) 

2/m-(i-l)(l)  ■ 

2/m-l(l) 

2/o(2) 

2/i(2) 

2/2(2) 

2/i.  (2) 

2/i  (2) 

2/2(2) 

2/3(2) 

2/i.+i(2) 

2/m-(i.  +  l  )(2) 

2/m  —  i,  (2) 

2/m-(L- ]  j(2) 

2/m—  1  (2) 

yo{B)  yi{N)  y2(-W) 

2/i i-P)  '  yi(N)  y3(N) 


Vl(N) 

Vl+i(N) 


[  y,n~(L+i){N)  ym-L(N )  ym-(L-i)(P)  •••  J 

Here  L  is  the  order  of  prediction  and  b  is  the  coefficient  vector  of  the  polynomial  B(z ) 
given  by 

B{z)  =  1  +  biz1  +  2  +  •  •  •  +  bizL .  (8) 

for  the  noiseless  case,  L  can  be  any  integer  greater  than  or  equal  to  the  model  order  n; 
however,  in  the  presence  of  noise  choosing  L  >  n  results  in  more  accurate  parameter 
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estimates  (see  Section  II).  Note  that  all  of  the  N  snapshots  are  used  simultaneously 
Lu  estimate  a  single  set  of  prediction  coefficients  (and  therefore,  a  single  set  of  poles). 

The  TLS-Prony  method  considers  the  effect  of  noise  perturbation  of  both  Y  and 
■y,  and  the  TLS  solution  attempts  to  minimize  the  effect  of  these  perturbations  on 
the  prediction  coefficient  vector  b  (see  [9,  10]  for  details).  This  is  accomplished  by 
obtaining  an  SVD  of  the  matrix  y  :  Y  and  truncating  all  but  the  first  n  sin¬ 
gular  values  to  arrive  at  an  estimate  y  ;  Y  [9,  10].  Inserting  y  :  Y  in 
Equation  (5)  gives  the  modified  linear  prediction  equation 

Yb=  -y  (9) 

from  which  the  linear  prediction  coefficient  vector  estimate  b  is  found  as 

b=-Y+y,  (10) 

where  +  denotes  the  Moore-Penrose  pseudoinverse.  A  numerically  robust  solution  for 
b  can  be  found  directly  from  the  SVD  of  y  ;  widehatY  >  as  is  shown  in  [10]. 
Finally,  the  estimates  for  the  poles  are  found  by 

Pj  =  zeroj  ,  j  =  1, 2, . . . ,  L.  (11) 

Once  the  L  poles  are  determined  from  Equation  (11),  one  must  separate  the  n 
“true”  poles  from  the  remaining  L  —  n  “extraneous”  or  “noise”  poles.  A  popular 
approach  is  to  choose  n  poles  based  on  their  location  with  respect  to  the  unit  circle. 
For  example,  one  can  choose  the  n  poles  closest  to  the  unit  circle  if  it  is  known 
that  the  poles  are  undamped  [29]  or  the  n  poles  with  smallest  moduli  if  it  is  known 
that  the  poles  are  damped  [•!].  However,  these  methods  do  not  apply  when  the  true 
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poles  may  lie  both  inside  and  outside  the  unit  circle.  In  this  case  we  can  classify 
poles  as  true  or  extraneous  based  on  the  energy  of  the  corresponding  mode.  We 
have  found  this  method  to  be  more  reliable  than  other  procedures  for  the  case  when 
the  true  modes  lie  both  inside  and  outside  the  unit  circle.  This  arises,  for  example, 
in  the  radar  scattering  problem  where  measurements  are  made  over  a  small  relative 
bandwidth,  and  the  exponential  modes  in  the  data  can  be  decaying  or  growing  over 
that  band  [3U,  31]. 

In  this  energy  criterion  method,  the  L  sets  of  amplitude  coefficients  can  be  found 
using  the  pole  estimates,  and  Equation  (3)  leads  to  the  following  least  squares  equa¬ 
tion  for  the  amplitude  coefficients, 


1  1 

pi  h 

Pi  Pit 

1 

Pl 

■■  Pl 

’  *id) 
i--'(l) 

*i(2)  ■ 

h(2)  ■ 

••  2i(Ar)  ’ 
■■  x-j(N) 

ss 

3/(1)  3/(2)  ■ 

i y(N) 

1  g"-1  • 

■  p r1 

.  *i( 1) 

2 l(2)  ■ 

■  ■  xL(N)  _ 

(12) 


or 

AlXl  »  Ya.  (13) 

A  least  squares  solution  to  Equation  (13)  can  be  computed  as 

XL  =  (AiALy'  AiY.  =  A*Y.<  (14) 


where  *  denotes  complex  conjugate  transpose.  (However,  in  practice,  more  numeri¬ 
cally  robust  procedures,  such  as  a  QR  decomposition,  should  be  used  to  solve  Equa¬ 


tion  (13).)  Because  only  n  singular  values  of 


V  ■  y 


are  nonzero,  there  are  at 
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most  n  pole  estimates  which  can  correspond  to  data  modes.  Therefore,  only  the  n 
poles  which  have  the  largest  energy  are  retained.  This  is  done  by  computing  the  L 
mode  energies  as 

N  m— 1 

E,  =  £  |Si(i)|2  £  N2’  i  =  1,2, ...,i  (15) 

t= 1  9=0 

and  retaining  those  n  poles  whose  corresponding  energies  are  highest.  We  then  rees- 
timate  the  amplitude  coefficients  of  these  n  poles.  This  is  done  using  Equation  (14), 
except  that  Ai  is  replaced  by  A ,  where  A  is  the  Vandermonde  matrix  composed  only 
of  the  n  columns  of  Ai  corresponding  to  the  n  selected  poles.  We  note  that  the  above 
procedure  produces  consistent  estimates  as  a  — >  0,  as  is  shown  in  the  Appendix. 

III.  Statistical  Analysis 

In  this  section  we  present  a  first  order  statistics  of  the  complete  set  of  parame¬ 
ter  estimates  obtained  in  the  TLS-Prony  method.  Assume  the  data  is  given  as  in 
Equation  (2).  Let  u;,  and  a,  be  the  angle  and  magnitude,  respectively,  of  each  pole 
/),.  thus  //,  -  .Similarly  let  q(f)  and  ;J(t)  be  the  angle  and  magnitude  vectors, 

respectively,  of  each  vector  of  amplitude  coefficients  x(t), 

1  r 

7(0  =  7i(0  72(0  •••  7n(0 

1  T 

m  =  [/mo  a(o  •••  &(oj  ’  (16) 

where 

T 

x{i)  =  /A(f)e^b)  ...  [in{t)e3^  •  (!') 
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Define  following  parameter  vectors: 


ex  = 

Bv  = 
0  = 


7T(1)  /?T(1)  7T(2)  /?T(2) 


•  •  •  uJn  O']  Ck2 

e 


r 


7t(W) 


(18) 


The  following  theorem  gives  the  first  order  approximation  of  the  probability  density 
function  (pdf)  of  0. 

Theorem  1:  Let  6  denote  the  TLS-Prony  estimate  of  9  which  is  given  by  the  n 
highest  energy  mode  estimates  found  in  Equations  (11)  and  (14).  Then  to  a  first  order 
approximation  (as  o  — »  0),  the  pdf  of  9  is  given  by  the  unbiased  Normal  distribution 


6  ~  M  (0 ,£,), 


(19) 


where 


u  — 


(7(1.1) 

f'(l.l)7';J-‘(l) 

V(\ .N) 

(/(1,A/)T^-J(.V) 

V(l) 

V(1  )Ta~l 

■  -T^(\)V(\,N) 

-T-bijvqi) 

<(A,1) 

u(N,\yr:r'{\) 

T(A:,A') 

TfA'.AO'VbA/) 

V(N) 

VINjTa-1 

■  -T-\N)U(N,N) 

7’-1(A')T(A’,A')7>-1(Af) 

-T-'(N)V(N) 

7'-1(A)V7(A')7q-1 

V*(l) 

V(l)7>-‘(1) 

V*(N) 

2 

ZTq-1 

-T-lV^(  1) 

T-'WwTe-^i)  ■ 

■  -T-'V^iN) 

7’a-1V^(A')7„-1(A) 

- T-'zZ 

T-'ZTa-1 

(20) 

where  •  and  •  in  Equation  (20)  are  real  and  imaginary  part  operators,  respectively, 
and  where  U(t,r ),  V(t),  and  Z  are  n  x  n  complex  matrices  which  depend  on  6 ,  L , 
and  in.  The  specific  iormulas  lor  these  entries  can  be  found  in  the  Appendix.  Tp(t) 
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and  Ta  are  diagonal  matrices  given  by 

Te(t)  =  dlas  (ft(0’ AW’  ’’ /Mo) 

Ta  =  diag  •  (21) 

\ai  «2  o.nJ 

Proof:  See  the  Appendix.  D 

Several  properties  of  the  covariance  can  be  derived  from  the  structure  of  E#.  Some 
of  these  properties  are  presented  in  the  following  corollaries. 

Corollary  1:  From  Eg  in  Equation  (20),  Cov  J  =  0,  and  Cov  (a>;,  c^)  = 

0. 

Proof:  Consider  the  blocks  of  Eg  containing  the  covariances  of  interest,  which  are 
given  by  U{t,  t )T$~'{t)  and  ZTa~1 .  From  Equation  (66)  in  the  Appendix,  it  can  be 
seen  that  U{t,t)  and  Z  are  Hermitian.  If  follows  that  the  diagonal  elements  of  U(t,t) 
and  Z  are  zero.  Since  Ty{t)  and  la  are  real,  diagonal  matrices,  the  diagonal  elements 
of  U(t,t)Tp~l(t)  and  ZTa _1  are  also  zero,  which  gives  the  desired  result.  □ 

Note  that  when  t  r,  U(t,r)  is  not  Hermitian,  so  the  diagonal  elements  of  U(t,r) 
are  not  zero.  Thus  it  is  not  in  general  the  case  that  the  magnitude  of  a \(t)  and  the 
angle  of  Xj(t)  are  uncorrelated  for  i  ^  j. 

Note  that  from  Equation  (20)  the  angle  variances  are  equal  to  the  magnitude 
variances  except  for  the  transformation  matrices  Tp{t)  and  Ta.  These  transformation 
matrices  can  be  eliminated  by  rescaling  some  of  the  parameters  in  6.  The  required 
rescaling  is  obtained  by  using  the  relative  magnitudes  of  the  poles  and  amplitude 
coefficients  as  the  estimated  parameters  instead  of  their  absolute  magnitudes.  That 
is,  define  the  estimate  to  be  as  in  Equation  (18),  but  with  5,  and  /3i(t)  replaced  by 
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the  relative  magnitudes  ^  and  We  then  note  that  the  Jacobian  transformation 
from  6  to  6\  is  given  by 

J  =  diag  (4, 7/3(1),  In,  7>(2), . . . ,  In,  Tp(N),  In,  Ta) .  (22) 

Corollary  2:  E#,  —  cov  (ffi)  is  given  by  Equation  (20)  with  all  Tp(t)  and  Ta 
matrices  replaced  by  identity  matrices.  It  follows  that  the  covariances  of  parameter 
angles  are  equal  to  the  covariances  of  the  corresponding  relative  magnitudes. 

Proof:  Immediate  from  the  fact  that  E#,  =  JEgJ  with  J  defined  in  Equation  (22). 

□ 

We  can  also  consider  a  reparameterization  of  6  in  which  real  and  imaginary  parts 
of  the  amplitude  coefficients  and  poles  are  considered  as  parameters.  Let  us  denote 
such  a  reparameterization  as  #2,  with  corresponding  covariance  matrix  which  would 
give  E*2 . 

Corollary  3:  Let  v  denote  a  complex  parameter,  which  is  either  a  pole  p,  or  an 
amplitude  coefficient  X{ (t).  Then  var(Re{i/})  =  var(Im{i^})  and  Re{i/}  is  uncorrelated 
with  lm{t/|. 

Proof:  The  result  can  be  obtained  by  applying  the  Jacobian  variable  transforma¬ 
tion  from  polar  to  rectangular  coordinates  to  E^.  This  transformation  is  straightfor¬ 
ward,  but  tedious,  and  not  presented  here.  □ 

Corollary  f:  E g  is  independent  of  the  absolute  phase  references  of  the  ampli¬ 
tude  coefficients  within  each  snapshot,  <p(t),  and  independent  of  the  absolute  phase 
reference  of  the  poles,  4>. 

Proof:  The  result  follows  by  examining  the  expressions  for  Tp(t)  and  Ta  in  Equa¬ 
tion  (21),  U(t,r),  V(t),  and  Z  in  the  Appendix,  and  noting  that  they  remain  un- 
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changed  if  Tx(t)  is  replaced  by  e^^Tx(t)  and  pi  is  replaced  by 


□ 


IV.  Examples 

In  this  section  we  present  a  set  of  examples  which  illustrate  the  performance  of  the 
TLS-Prony  method.  We  first  compare  the  first  order  statistics  presented  above  to  the 
CRB  for  a  number  of  cases.  The  CRB  for  this  data  model  is  presented  in  [32].  We  then 
compare  the  first  order  statistics  to  those  obtained  using  Monte-Carlo  simulation. 

A.  Exam  pit  1:  Single  Exponential  Mode 

In  this  example  we  consider  a  single  pole  model  with  one  snapshot  of  data  (and 
thus  one  amplitude  coefficient).  The  experiment  entails  moving  the  pole  along  the 
positive  real  axis  from  0.1  to  10,  he.  0.1  <  p  <  10  (the  results  are  independent  of  the 
pole  angle  by  Corollary  4,  so  an  angle  of  zero  is  chosen).  For  each  pole  location  ,  W'e 
calculate  the  parameter  variances  using  Equation  (20)  for  data  sets  of  lengths  2,  5, 
10,  20,  50,  and  100.  For  comparative  purposes,  the  amplitude  coefficient  associated 
with  the  pole  is  chosen  to  be  a  positive  real  number  such  that  the  mode  energy 
( x 2  XTfco1  P 21 )  is  un>ty  for  each  pole  location  and  data  length.  The  noise  power  is  also 
kept  constant  at  a  —  1.  The  model  order  L  is  chosen  to  be  one  third  of  the  data 
lenglli  in.  which  has  been  shown  to  be  near  optimal  for  a  number  of  cases  (see  [18,  33], 
and  ^4.3  below). 
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1.  Pole  Variances 


The  first  order  theoretical  variances  for  the  estimated  pole  angle  and  pole  magni¬ 
tude  appear  as  the  dashed  lines  in  Figures  1  and  2,  respectively;  the  corresponding 
CRBs  appear  as  the  solid  lines  in  these  figures.  From  Figure  1  we  see  that  the  pole 
angle  variances  are  close  to  the  CRBs  when  the  pole  is  inside  the  unit  circle.  When 
the  pole  is  outside  the  unit  circle,  the  variances  become  much  higher  than  the  CRBs 
(except  for  the  rn  =  2  case).  For  larger  data  lengths  the  disparity  with  the  CRB  is 
much  more  pronounced.  This  is  because  backward  linear  prediction  is  used  in  our 
TLS-Prony  estimation  method.  With  backward  linear  prediction,  extraneous  poles 
lie  outside  the  unit  circle,  thus  making  estimation  of  poles  outside  the  unit  circle  more 
difficult  [34].  The  use  of  forward  linear  prediction  would  give  corresponding  results 
for  poles  inside  the  unit  circle.  Similar  observations  apply  to  the  pole  magnitude 
variances  (see  Figure  2).  The  pole  magnitude  variances  can  be  normalized  to  give 
relative  error  of  the  pole  magnitude,  i.e.  var^^-j.  If  this  is  done,  one  obtains  exactly 
the  same  curves  as  in  Figure  1  (c/.  Corollary  2). 

From  these  two  figures  we  also  see  that  inside  the  unit  circle  the  variances  for  pole 
angle  are  higher  than  the  variances  for  pole  magnitude  and  vice-versa  outside  the 
unit  circle.  This  is  because  the  angular  uncertainty  becomes  greater  as  a  pole  moves 
closer  to  the  origin.  As  expected,  the  pole  angle  variance  approaches  infinity  as  the 
pole  approaches  the  origin. 

From  Figures  1  and  2  we  see  that  the  pole  angle  and  magnitude  variances  are 
asymptotically  (as  m  — >  oo)  lowest  when  the  pole  is  on  the  unit  circle,  and  that  on 
the  unit  circle  the  variances  are  decreasing  by  1/m2  (m  is  the  data  length).  This  is 


consistent  with  the  well-known  1  /m3  variance  decrease,  since  the  amplitude  coefficient 
is  adjusted  in  this  experiment  to  keep  the  mode  energy  constant  (if  the  amplitude 
coefficient  is  left  unchanged,  the  variance  decrease  is  1/m3). 

When  the  pole  is  not  on  the  unit  circle,  the  variances  do  not  decrease  to  zero 
as  m  — >  oo.  Recall  that  we  keep  the  total  mode  energy  constant.  For  a  decaying 
or  growing  exponential  mode,  adding  data  points  while  keeping  the  energy  constant 
results  in  adding  data  points  with  smaller  and  smaller  amplitude.  As  a  result,  the 
parameter  estimate  variances  do  not  continue  to  decrease. 

2.  Amplitude  Coefficient  Variances 

1  he  variances  lor  the  amplitude  coefficient  angle  and  magnitude  appear  in  Fig¬ 
ures  3  and  4.  As  before,  each  curve  is  a  plot  of  variance  versus  pole  radius  for  various 
data  lengths  m.  There  are  several  points  to  note  in  these  figures.  First,  when  the  pole 
is  inside  the  unit  circle,  increasing  the  number  of  data  points  provides  no  significant 
decrease  in  the  variances.  The  first  data  point  with  no  noise  added  is  precisely  the 
amplitude  coefficient.  When  the  pole  is  inside  the  unit  circle,  this  amplitude  does 
not  change  much  as  a  function  of  data  length,  and  consequently  its  variance  does 
not  change  much.  However,  when  the  pole  is  outside  the  unit  circle,  the  amplitude 
coefficient  decreases  rapidly  toward  zero  as  data  length  increases.  Thus,  outside  the 
unit  circle  the  estimate  of  the  amplitude  cannot  be  expected  to  vary  much  around 
zero  and  the  magnitude  variances  become  low.  Also,  variance  of  the  estimated  angle 
becomes  quite  large  because  of  high  angular  uncertainty  for  points  near  zero. 
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3.  Prediction  Order  Considerations 


In  the  above  figures  we  used  a  prediction  order  L  equal  to  one-third  the  data 
length.  We  next  consider  the  effect  of  prediction  order  on  the  variances  of  the  TLS- 
Prony  parameter  estimates.  We  consider  a  single  exponential  mode  whose  pole  is  on 
the  unit  circle.  We  choose  a  =  1  and  choose  the  amplitude  coefficient  so  that  the 
mode  energy  is  unity,  as  before. 

Figures  5,  6,  and  7  show  the  variances  of  the  pole  angle,  amplitude  coefficient  mag¬ 
nitude,  and  amplitude  coefficient  angle.  The  solid  lines  represent  the  CRBs  (which 
are,  of  course,  independent  of  TLS  prediction  order),  and  the  dotted  lines  represent 
the  TLS-Prony  variances.  Figure  5  has  been  presented  in  earlier  works  [18,  23,  25], 
but  the  amplitude  coefficient  was  not  considered  there.  Since  the  pole  is  on  the  unit 
circle,  the  pole  magnitude  results  are  identical  to  the  pole  angle  results  (c/.  Corol¬ 
lary  2,  with  q j  =  1).  From  these  figures  we  can  see  that  the  best  prediction  order 
choice  is  approximately  rn/3]  this  agrees  with  the  results  in  [18,  23,  25]. 

B.  Example  2:  Two  Undamped  Exponential  Modes 

In  this  example  we  consider  two  poles  at  a1e-du'0+Au'/2)  and  c^e-7^0-^/2)  with 
O]  =  a2  =  1.  Variances  are  computed  for  various  data  lengths  (m  =  5,  10,  20,  50, 
and  100)  and  angle  separations  Aw.  The  variance  results  are  independent  of  uo  so 
uq  =  0  is  used.  Again,  L  =  m/3,  a  =  1,  one  snapshot  of  data  is  used  ( N  =  1),  and 
each  amplitude  is  chosen  such  that  the  corresponding  mode  energy  is  unity. 

Figure  8  shows  the  variances  for  the  pole  angle  estimate  (of  either  pole)  as  a 
function  ol  pole  separation  and  data  length.  The  solid  lines  again  show  corresponding 
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CRB  results.  The  variances  for  the  pole  magnitudes  are  equal  to  the  pole  angle 
variances  because  these  poles  are  located  on  the  unit  circle.  We  see  that  the  TLS- 
Prony  algorithm  performs  well  with  respect  to  the  CRB.  In  fact,  for  the  curves  shown, 
the  TLS-Prony  variance  curves  are  always  within  2dB  of  the  corresponding  CRB 
curves. 


C.  Example  3:  Monte-Carlo  Simulation  of  Ten  Mode  Case 

In  this  example,  we  have  chosen  ten  exponential  modes  to  represent  a  general 
case.  The  true  pole  location  of  each  mode  is  denoted  with  an  “x”  in  Figure  9(a). 
For  this  case  we  have  rn  =  100  data  points,  L  =  33,  and  a  =  0.001.  The  amplitude 
coefficients  are  chosen  so  that  each  mode  energy  is  unity;  this  corresponds  to  an  SNR 
of  lOdB  per  mode;  the  total  SNR  (signal  power/noise  power)  is  20dB.  The  modes 
were  chosen  as  a  representative  sampling  of  various  situations. 

Figure  9  presents  a  comparison  between  the  TLS-Prony  estimate  theoretical  vari¬ 
ance's  and  variance  estimates  obtained  using  Monte-C’arlo  simulations.  The  theoreti¬ 
cal  variances  are  shown  as  two-standard  deviation  concentration  ellipses  around  each 
pole.  These  ellipses  (they  are  actually  circles,  by  Corollary  3)  give  87%  confidence 
intervals  in  the  complex  plane  for  pole  pole  estimates.  The  dots  in  Figure  9  are  pole 
estimates  corresponding  to  the  ten  highest  energy  poles  from  each  of  100  Monte-Carlo 
simulations.  The  details  of  the  pole  estimation  accuracy  are  summarized  in  Table  1. 
Note  that  the  Monte-Carlo  variances  are  close  to  those  predicted  by  theory  for  most  of 
the  poles,  in  particular  for  those  closer  to  the  unit  circle.  For  poles  well  inside  the  unit 
circle,  there  is  some  bias  present  which  is  not  predicted  by  a  first  order  approxima¬ 
tion;  in  addition,  the  predicted  variance  is  smaller  than  the  actual  variance.  As  poles 
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Table  1:  Theoretical  and  simulation  variances  for  the  poles  of  a  general  ten  mode 
case.  ■ _ 


Pole 

number 

Theoretical  variance 
(dB) 

Simulation  variance 
(dB) 

Pole 

radius 

1 

-24.2 

-22.4 

1.115 

2 

-42.7 

-43.4 

1.05  ; 

3 

-35.7 

-35.2 

0.8 

4 

-25.8 

_ 

-22.4 

1.12 

5 

-57.1 

-56.9 

1.0 

6 

-57.0 

-56.1 

1.0 

7 

-25.8 

-21.5 

0.6 

8 

-29.4 

-29.1 

1.115 

9 

-39.1 

-32.2 

0.9 

10 

-30.1 

-26.7 

0.7 
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move  outside  the  unit  circle  to  the  radius  of  the  extraneous  poles,  some  deterioration 
occurs  in  terms  of  misidentifying  pole  estimates  as  “true”  or  “extraneous”.  Note  the 
rapid  increase  in  variance  of  a  pole  estimate  as  its  radius  increases,  by  comparing  the 
variance  for  pole  2,  8,  4,  and  1. 


V.  Conclusions 

We  have  presented  a  statistical  analysis  for  estimated  poles  of  the  TLS-Prony 
algorithm.  This  analysis  includes  complete  expressions  for  the  covariance  matrix  of 
the  parameters  of  an  exponential  model  which  contains  one  set  of  poles  and  multiple 
sets  of  amplitude  coefficients.  The  poles  of  this  model  may  lie  anywhere  in  the  complex 
plane.  Using  these  expressions,  several  useful  properties  of  the  covariance  matrix  were 
established.  These  include  independence  of  the  two  parameters  for  each  amplitude 
coefficient  and  pole,  whether  one  considers  a  polar,  a  relative  magnitude  polar,  or  a 
rectangular  real  and  imaginary  part  parameterization.  It  was  also  established  that 
the  variances  of  these  pairs  are  equal  for  the  relative  magnitude  polar  and  rectangular 
real  and  imaginary  part  parameterizations. 

The  results  of  this  paper  provide  a  sound  basis  for  performance  analyses  of  the 
TLS-Prony  estimation  method.  We  have  extended  previous  works  to  include  the 
general  damped  undamped  case,  as  well  as  to  include  amplitude  coefficient  parameters 
in  addition  to  pole  parameters.  The  results  can  be  used  to  analyze  various  situations 
and  evaluate  the  potential  success  of  applying  the  TLS-Prony  estimation  algorithm, 
as  the  corollaries  and  examples  in  the  paper  demonstrate. 
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Appendix:  Proof  of  Theorem  1 

From  Equation  (10)  we  note  that 

6  =  fc_6  =  - Y+y  +  S+s 

=  -Y+(s  +  s)  +  S+s 

=  ~(y+~  S+)  s  -  Y+s  (23) 

where 

y  :  Y  —  s  :  S  +  s  :  S  ■  (24) 

Here,  s  :  5  is  the  noise  free  version  of  y  :  y  ,  and  6  is  the  Lth  order  linear 
prediction  coefficient  vector  associated  with  the  poles  of  the  noiseless  model.  Note 
that  6  gives  the  n  true  poles  and  L  —  n  extraneous  poles  (c/.  [34]).  In  order  to  evaluate 
the  expression  in  Equation  (23)  we  use  the  following  identity  for  any  matrix  M  [35] 

M+  -  M+  =  -M+MM+  +  (M*M)+  M'P^  +  Pi.AT  (MM')+  ,  (25) 

where  M  =  M  -  M ,  =  /m  -  MA/+,  P4.  =  /„  -  +  ,  and  lq  is  the  q  x  9 

identity  matrix.  Using  the  fact  that  PgS  —  0  and  a  first  order  approximation,  we 
then  obtain 

6  -  -Y+Sb-  P+tS~S'+S+s-Y+s 

«  -5+56  -  Ps.S*S*+S+s  -  S+s.  (26) 
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The  above  approximation  is  valid  since  the  matrices  y  :  y  and  $  :  S  have 


the  same  rank.  Let 


y  :  Y  —  s  :  S  +  w  :  W  > 


where  w  :  W  is  defined  as  the  noise-only  part  of  y  :  Y  (see  Equations  (2) 
and  (7)).  By  using  the  perturbation  analysis  in  [36]  on  the  matrices  s  :  S  ■> 
y  .  y  ,  and  y  :  y  ,  it  can  be  shown  that  to  a  first  order  approximation 

S+  5:5  =*S’+  w  :  W  ■  (28) 


The  above  equation  also  implies  that  S"Sm+  =  VE*5*+.  Thus,  Equation  (26)  can  be 
written  (to  a  first  order  approximation)  as 

b  =  -S+Wb  -  Pg.W*S*+ S+  s  -  S+w.  (29) 


From  Equation  (29),  we  note  that  b  -*  0  as  a  -*  0  since  the  elements  of  W  and  w  are 
uncorrelated,  zero  mean,  complex  white  Gaussian  random  variables.  Therefore,  the 
resulting  L  pole  estimates  are  consistent  as  cr  — >  0.  Similarly  it  can  be  shown  that 
the  L  sets  of  the  amplitude  coefficients  are  also  consistent  as  a  — *  0.  Note  that  the 
“true”  amplitude  coefficients  of  the  extraneous  modes  are  zero;  thus  it  follows  that 
choosing  the  n  highest  energy  poles  as  the  true  poles  is  consistent  as  cr  — >  0. 

Note  that  SR j.  =  0.  Multiplying  both  sides  of  Equation  (29)  through  by  5,  we 
obtain 

Sb  =  -55+e,  (30) 

where  t  —  w  +  Wb. 
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From  Equation  (2),  we  can  write  S  as 


-ti  (l)Pi 


*2(1>2 


( t  )Pn 


*i(1  )P?~(L  +  1)  *2(1  )p?~(L+1) 


*l(2) 

*i(2)pi 


®2(2)P2 


•••  *n(l)pr(i+1) 


^n(2)pn 


x1(2)p”-(i+1)  r2(2)p^-(L+1) 


x„(2  )p^~(L+1) 


*i{N) 

x1{N)pl 


x2(N) 

x2(N)p2 


xn(N) 
Xn  (^)Pn 


*iW 


*2(iV)p? 


xn(N)p„ 


m-(L+ 1) 


where  G  is  given  by 


5  =  HG , 


Equation  (30)  thus  becomes 


Pi  Pi  •••  Pi 


P2  p\  •••  P2 


Pn  Pn  •••  Pn 


/7G6  -  -HGS+t. 


Now  note  from  Equation  (8)  that  the  true  and  estimated  Lth  order  characteristic 


polynomials  are  B{z)  =  l+61z'+62z2+ - \-bLzL  and  B(z)  =  l+biz1+b2z2-i - \-bLzL, 


respectively.  Note  that  B  ( pi )  —  0  and  B  (p,)  =  0. 

We  can  use  a  first  order  Taylor  expansion  to  find  an  expression  for  the  error  in 
the  estimated  pole  locations.  We  follow  the  technique  in  [17].  For  each  pi  we  obtain 


0  =  B{pt) 

=  B  (pt)  +  ~B(z)\i=Pi  (pi  -  p,)  +  (higher  order  terms) 

=  B  (pi)  -  B  (pt)  +  —  £(2)|2=p,  (P‘  _  P*')  +  (higher  order  terms) 

~  1  +  Wp,  +  h P2t  +  ■■■  +  l>LPi  -  (l  +  biPi  +  b2p 2  + - f  bLpfj 


Pi  Pi  ■■■  Pi 


+  LlLpf-') 

(Px  ~  Pi) 

bi  -  bi 

62  —  62 

1 - 1 

-O 

CM 

-O 

-O 

+ 

1 

2 Pi 

bL  ~  bL 

LPi  M  _ 

(b-b)+  Vi  (pi  -  pi) . 

Thus,  to  a  first-order  approximation 

(Pi-Pi)  =  ~-  pi  p\  ■ 
1 h 

where  pt  is  given  by 


•  pf  b , 


(Pi -Pi) 


7h  -  bi  b2  ■  ■  ■  bL 
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In  matrix  form  we  thus  obtain 


P=P-P=  -FGb, 


(38) 


where 


P  = 
P  = 


iT 


Pi  P2  ■■■  Pn 


Pi  P2  Pn 


F  =  diagfil, 

\Vl  V2  VnJ 


(39) 


Since  the  n  poles  are  distinct,  H  is  of  full  column  rank.  Hence,  we  can  multiply 
Equation  (34)  by  ( H’H)~ 1  H’  to  get 


Gb  -  -GS+t, 


and  by  substituting  Equation  (40)  into  Equation  (38)  we  obtain 


(40) 


P  =  FGS+t.  (41) 

We  now  note  that  to  a  first  order  approximation,  P  is  given  by 

P  =  a  ©  eJW  —  a  0  eJUJ 

=  (a  +  5)0ej(w+;)-o0eJ“ 

=  (a  -f  5)  ©  eJU/  ©  (1  +  jui  +  (h.o.t.))  —  a  Q  eJUJ 

»  ( Q  +  Q )  ©  cjw  +  Q  ©  e-JW  ©  ju)  —  q  ©  t3UJ 

=  Tp1  {Ta5  +  ju>),  (42) 
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where  0  denotes  the  Hadamard  product,  and 

r  -.7 

eJW  =  gi^l  .  .  .  eju> n 

1T 

u  —  U>1  U>2  •  ■  •  U>n 

T 

a  ~  «i  a2  •••  5n 

Tp  =  diag  )  (43) 

\Pl  P2  PnJ  ’ 

and  where  u,  —  £>,•  —  ut  and  a,  =  a,-  —  a,-.  From  Equations  (41)  and  (42)  we  obtain: 

w  =  lm{TpFGS+e} 

5  =  Re  {T“iTpFG’5'+f}  .  (44) 

Recall  that  the  elements  of  W  and  w  are  uncorrelated,  zero  mean,  complex  white 
Gaussian  random  variables.  Thus,  t  is  multivariate  Gaussian  with  zero  mean  and 
covariance  matrix 

E[t(']  =  E  (  w  ;  W  *  )  (  w  :  W  =oDD\  (45) 

D  is  defined  as  a  (in  —  L)  N  x  mN  block  diagonal  matrix  given  by 

D  =  diag  (B,  B,...,  B) ,  (46) 
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where  B  is  given  by 


1  bx  b2  ■  ■  ■  bL  0  0  •  •  •  0 

0  1  &!  •••  bL  0  •••  0 

;  :  .  (47) 

0  0  1  6,  •••  i>i_,  bL  0 

0  •••  0  0  1  •••  bi- 2  6l-  i  bi 

-*  (m— L)x(m) 

We  also  have 

e[ut]=E  ([„,  ;  W  ][[])([»  :  w][j])  =0.  (48) 

Using  these  results,  along  with  the  following  relationships  (proven  in  [16]) 

Re{u}Re{v7}  =  -  |Re{unr}  -f  Re{nt>*}j 
Im{u}lm{t>7}  =  — -  [Re{uuT}  —  Re{un’')] 

Re{u}lm{u7  }  =  -  [lm{unr}  —  lm{nn*}  ,  (49) 

we  obtain  the  following  covariances  for  the  pole  parameters: 

Epu7}  =  ^R  e{TpFGS+DD'S'+G'F*T;} 

Fpa7}  =  {tvfgs+dd'S’+g*f't;t--1} 

F  [dm7]  =  -Re{T:%FGS+DD'S'+G’F'T;T'a-').  (50) 

To  obtain  the  covariances  for  the  amplitude  coefficient  parameters  we  use  Equa¬ 
tion  (14),  which  provides  the  amplitude  coefficient  estimates  for  each  snapshot  in 
terms  of  the  estimated  poles.  We  now  note  the  following 
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(51) 


=  A+{Sa  +  Na)-A+Sa 
=  (A+  -A+)Sa  +  A+Na, 

where 

x=  1(1)  x(2)  •••  x(N)  ,  (52) 

and  where  Sa  is  the  noise  free  version  of  Ya  and  Na  =  e(l)  e(2)  •••  e(iV)  is 
the  corresponding  noise. 

We  apply  the  identity  in  Equation  (25)  to  the  first  term  of  Equation  (51).  Since 
the  n  poles  are  distinct,  A  and  A  are  full  rank.  Also,  since  m  >  n  and  Sa  6  Range(A), 
we  have 

(A+  -  A+)Sa  =  -A+AA+Sa.  (53) 

From  Equations  (51)  and  (53)  we  obtain  the  following  first  order  approximation 

X  =  -AJrAA+Sa  +  A+Na 
=  -A+AX  +  A+Na 

«  -A+AX  +  A+Na.  (54) 

Note  that  A  is  a  matrix  in  which  each  column  is  composed  of  the  amplitude  coefficient 
variations  for  each  snapshot: 

A' =  £•(])  x{2)  •••  x(N)  ■  (55) 

Following  the  same  procedure  as  in  Equation  (42),  a  first  order  approximation  of  x(t) 
is  given  by 

(56) 
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where 


T 

7(0  =  7l(0  72(0  •••  7n  (0 

~  r  T 

m  =  [aw  a(<)  •••  &(t)_ 

T*(i)  =  diag(j7w’^’"-’^)  (57) 

and  where  7,(f)  =  7 ,-(<)  -  7,(f)  and  /?,(*)  =  #(*)  -  #(*)■  From  Equations  (54),  (55), 
and  (56)  we  obtain 

7(0  =  lm{Tx(t)A+  (-Ax(t)  +  e(t))| 

P(t)  =  R  e[Tp\t)Tx{t)A+ (-Ax(t)  + e(t))}  .  (58) 

Beiore  computing  covariances  for  the  parameters  in  Equation  (58),  we  need  to 
perform  some  manipulations  to  Ax(t),  since  the  random  variable  A  does  not  appear 
at  the  rightmost  position.  We  proceed  with  a  first  order  approximation  as  follows 


Ax(t)  — 


Pi 

-  pi 

P2 

~P2 

Pn 

~Pn 

Pi 

-pi 

P\ 

-Pi 

••  Pi 

“Pn 

p'r 1 

-  p?-1 

pT1 

1 

“e 
m  3 

7 

fir1 

“Pn 

0 

0 

Pi 

P2 

‘2p2p2 


0 

Pn 

2  PnPn 


{rn  -  ] )  p\n  2pi  (m  -  1 )  p™  2p2  ■  ■■  (m  -  1 )  p”1  2pn 
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=  CATpdi&g  (pi,P2,  ■  ■  -iPn)  x{t) 


=  CATx'\t)TpP 

=  CAT~\t)TpFGS+t,  (59) 

where  C  is  a  diagonal  matrix  given  by 

C  =  diag  (0, 1, . . .  ,m  —  1) .  (60) 

Equation  (58)  can  now  be  approximated  by 

7(0  ss  hn{Tx(t)A+  (-CAT;\t)TpFGS+t  +  e(i))} 
m  ~  Re{T^(t)Tx(t)A+  (~CAT-1(t)TpFGS+e  +  e(t))}  ,  (61) 

Since 

E[te*(t)]  =  E^w  w  ]  [|j)e*w]  =*D(i),  (62) 

where  D(t)  are  each  given  by  the  tth  column  block  of  D  ( cf .  see  Equation  (46)),  we 
also  have 

E{tfT(t)\  =  £■[([„;  ^  ]  []J)eT(0]  =0 
E[e(t)er(i)}  =  0.  (63) 

Using  these  results  and  Equations  (49)  and  (61)  we  obtain  the  following  covari¬ 
ances  for  the  amplitude  coefficient  parameters 

E  [7(07T(r)]  =  yRe{Tx(t)A+  (C AT;\i)TPFGS+ DD* S+*G* F‘T;T;-' (r)A*C* 
- C ATpT~i{t)FGS+D(r)  -  D’{t)S+’GmF'T;T;-l{r)A-C* 

+  lm6t,r)A+'T;(r)} 
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E  [i(t)pT{r)}  =  ^lm{Tx(t)A+  (< C  AT;\t)TvFGS+ DD*  S+*G*  * 

—CAT~l(t)TpFGS+D(r)  -  D'{t)S+'G' F'T;T'x~\r)A'C * 

+Jm6t,r)  A+'T*{r)T^’{r)} 

F  \p{t)pr{r)}  =  ^Re{Tp\t)Tx{t)A+  [C  AT;\t)TpFGS+ DD' S+' G' F'T;T;-\r)A'C' 
-C AT~1(t)TpFGS+  D(r)  -  D'{t)S+'G'F'T;T*~\r)A'C' 
+Im6t,r)A+*T;(r)Tp-''(r)}.  (64) 

Using  Equations  (44),  (49),  and  (61)  we  can  also  compute  the  cross-covariances 
between  the  poles  and  the  amplitude  coefficients  as  follows 


E  [7(0^T] 
E  [7 {t)aT 

e  \m*T\ 

F  p(0or] 


|r e{Tx{t)A+  ( -CAT-\t)TpFGS+DD '  +  £>*(<))  S+’G'F'T;} 

|lm  {Tx(t)A+  ( -CATx-\t)TpFGS+DD *  +  £>*(<))  5+*G*FT;r;1‘} 

~ ^Im  ^T01(i)Tx(t)A+  ( -CAT;\i)TpFGS+DD '  +  £•(*))  5+’GTT;} 
^Re{77'(/)7;.(0/l+  {-CAT;\t)TpFGS+DD-  +  Dm(t))  S^G'F'T;T^'}. 

(65) 


Equations  (50),  (64),  and  (65)  completely  specify  £0  as  given  in  Equation  (20) 
using  the  following  substitutions 

U(t,r)  =  R(t)Z  R'(r)  -  R(t)Q(r)  -  Q'(t)R'(r)  +  ^Tx{t)A+ A+mT'x{r)8UT 
V(t)  =  -R(t)Z  +  Q'(t) 

Z  =  ^TPFGS+DD' S' + G'  F'T'p 
Q(t)  =  a-TpFGS+D{i)A+'T'x[t) 

R(t)  =  Tx(t)A+C  AT~\t).  (66) 
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Furthermore,  since  t  and  e(t)  are  zero  mean,  Equations  (44)  and  (61)  imply  that 
the  parameter  estimates  are  unbiased,  i.e.  to  a  first  order  approximation  E  6  =  6. 
This  completes  the  proof.  □ 
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pole  magnitude  variances 


Figure  2:  Pole  magnitude  variances  for  single  pole  data  (n  =  1). 
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amplitude  coefficient  angle  variances 


figure  3:  Amplitude  coefficient  angle  variances  for  single  pol 
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amplitude  coefficient  magnitude  variances 


Figure  4:  Amplitude  coefficient  magnitude  variances  for  single  pole  data  (n=l). 
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amplitude  coefficient  magnitude  variances 


prediction  order 


Figure  7:  Amplitude  coefficient  angle  variances  for  various  prediction  orders. 


Figure  9:  A  general  ten  mode  case:  (a) 
deviation  theoretical  bounding  circles  and 


B.  Preprint  of  “A  Modified  TLS-Prony  Method  Using  Data  Decimation 
for  Computationally  Efficient  Estimation  of  Damped  Exponentials  in 

Noise” 


The  following  pages  contain  a  preprint  of  the  paper,  “A  Modified  TLS-Prony  Method 
Using  Data  Decimation  for  Computationally  Efficient  Estimation  of  Damped  Ex¬ 
ponentials  in  Noise,”  which  has  been  submitted  to  IEEE  Transactions  on  Signal 
Processing. 


58 


A  Modified  TLS-Prony  Method  Using  Data  Decimation  t 


William  M.  Steedly  Ching-Hui  J.  Ying  and  Randolph  L.  Moses 

The  Analytic  Sciences  Corporation  Department  of  Electrical  Engineering 
Reston,  VA  22090  The  Ohio  State  University 

Columbus,  OH  43210 


Abstract 


This  paper  introduces  a  modified  TLS-Prony  method  which  incorporates  data  dec¬ 
imation.  The  use  of  data  decimation  results  in  the  reduction  in  the  computational 
complexity  because  one  high  order  estimation  is  replaced  by  several  low  order  estima¬ 
tions.  We  present  an  analysis  of  pole  variance  statistics  for  this  modified  TLS-Prony 
method.  This  analysis  provides  a  quantitative  comparison  of  the  parameter  estima¬ 
tion  accuracy  as  a  function  of  decimation  factors.  We  show  that  by  using  decimation, 
one  can  obtain  comparable  statistical  performance  results  at  a  fraction  of  the  com¬ 
putational  cost,  when  compared  to  the  conventional  TLS-Prony  algorithm. 
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I.  Introduction 


A  popular  high  resolution  estimation  technique  is  the  use  of  backward  linear  pre¬ 
diction  coupled  with  singular  value  decomposition  (SVD)  and  total  least  squares  [1], 
here  called  the  TLS-Prony  technique.  This  technique  has  been  shown  to  provide 
good  parameter  estimates  of  damped  exponential  signals  in  noise  for  various  types 
of  data  [1,  2].  However,  for  large  data  lengths,  the  TLS-Prony  method  can  be  com¬ 
putationally  expensive.  The  reason  for  this  is  that  the  TLS-Prony  method  involves 
computing  the  Singular  Value  Decomposition  (SVD)  of  a  data  matrix  of  size  (m,n) 
where  m  is  related  to  data  length  and  n  to  prediction  order.  For  best  accuracy  of  the 
parameter  estimates,  n  ~  y.  Thus,  this  data  matrix  becomes  quite  large  for  large 
data  lengths.  To  overcome  this  problem  it  is  sometimes  possible  to  decimate  the  data 
before  applying  the  TLS-Prony  technique;  the  result  is  often  a  large  reduction  in 
computations.  In  this  paper  we  consider  the  statistical  and  computational  properties 
of  the  TLS-Prony  algorithm  when  used  in  conjunction  with  data  decimation. 

Data  decimation  has  been  considered  before  [3,  4]  in  the  context  of  spectral  esti¬ 
mation.  This  technique  entails  using  only  part  of  the  measured  data.  Decimation  of 
correlation  sequences  was  also  considered  in  [5];  this  technique  effectively  uses  all  the 
measured  data,  but  is  somewhat  restrictive  in  that  it  applies  only  to  correlation-based 
parameter  estimation  techniques.  These  works  do  not  present  a  quantitative  analysis 
of  statistical  properties  of  the  resulting  parameter  estimates. 

In  this  paper  we  develop  a  data  decimation  technique  which  is  based  on  the  TLS- 
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Prony  algorithm  [1],  We  also  present  a  theoretical  statistical  analysis  of  the  accuracy 
of  the  TLS-Prony  parameter  estimates  when  decimation  and  any  linear  FIR  filter 
are  used.  Based  on  this  analysis,  we  present  a  quantitative  comparison  of  estimation 
accuracy  for  various  types  of  data  decimation  schemes.  In  particular,  we  compare 
the  decimated  and  nondecimated  procedures  in  terms  of  estimation  accuracy.  Our 
analysis  demonstrates  that  the  performance  of  the  decimated  TLS-Prony  procedure 
is  comparable  to  the  performance  of  the  nondecimated  TLS-Prony  procedure  for  un¬ 
damped  exponential  modes.  We  also  develop  a  complexity  analysis  and  show  that  the 
decimated  algorithm  is  computationally  more  efficient  than  the  nondecimated  algo¬ 
rithm.  Both  the  statistical  performance  and  the  decrease  in  computational  complexity 
are  verified  by  Monte-Carlo  simulations. 

We  then  apply  the  statistical  analysis  to  consider  two  specific  cases  of  interest. 
First,  the  signals  of  interest  may  be  bandlimited  and  occupy  a  relatively  small  region 
of  the  unambiguous  frequency  range  /  €  [—  § ,  §  -  In  this  case  one  is  interested  in 
analyzing  a  subset  of  the  whole  frequency  range.  For  example,  this  technique  was  used 
to  investigate  radar  signatures  of  aircraft  [2,  6].  One  question  is  whether  decimation 
can  improve  the  accuracy  of  the  pole  estimates  in  this  case. 

A  second  case  of  interest  is  when  the  signal  occupies  most  or  all  of  the  unambiguous 
frequency  range.  In  this  case  we  filter  the  data  to  isolate  a  number  of  subbands,  then 
use  a  decimation  version  of  TLS-Prony  to  estimate  the  poles  in  each  of  the  subbands. 
This  idea  is  similar  in  principle  to  beamspace  prefiltering  in  array  processing  [7,  8,  9]. 
By  focusing  on  particular  bands  one  at  a  time,  estimation  techniques  can  be  used 
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with  lower  model  orders  since  there  are  typically  fewer  modes  within  each  of  the 
bands.  Thus,  a  single  wideband  estimation  procedure  is  replaced  by  several  lower 
order  estimations.  This  has  the  advantages  of  being  much  less  numerically  intensive 
and  of  being  amenable  to  parallel  implementation  [9]. 

An  outline  of  this  paper  is  as  follows.  In  Section  II  we  develop  the  modified 
TLS-Prony  procedure.  In  Section  III  we  derive  the  first  order  approximation  of  the 
statistics  of  the  estimated  parameters.  In  Section  IV  we  develop  a  procedure  and  a 
complexity  analysis  for  performing  full  spectrum  estimation.  We  also  discuss  filter 
design  and  performance  loss  in  the  estimation.  Section  V  presents  some  simulation 
studies  using  decimation.  Finally,  Section  VI  concludes  the  paper. 

II.  Decimation  Estimation  Procedure 

A.  Data  Model 

Assume  we  have  N  “snapshots”  of  data  vectors  y(t),  each  of  length  m : 

1T 

y{i)=  yo(t)  yi{t)  ■■■  ym-i(t)  t  =  l,2,...,N.  (1) 

Each  data  vector  is  modeled  as  a  noisy  exponential  sequence 

n 

yq{i)  =  J2xi(t)Pi  +  q  =  0,1,...,  m  —  1.  (2) 

i=l 

There  are  n  exponential  modes  in  the  data;  the  n  poles  {p,}"=1  do  not  vary  from 
snapshot  to  snapshot,  but  the  amplitudes  Xi(t)  may  vary.  Here,  it  is  assumed  that 
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(eg(t)}  are  uncorrelated  zero  mean  complex  white  Gaussian  noise  sequences  with 
variance  a.  Equation  (2)  may  be  compactly  written  as 


y{t)  =  Ax(t)  +  e(t),  (3) 

where  e(t)  =  e0(t)  e^t)  •••  etn_1(<)  j  ,  *(*)  =  Xo(t)  Zl(t)  •••  xn^(t)  , 

and  A  is  the  m  x  n  Vandermonde  matrix  derived  from  n  signal  poles 

1  1  •••  1 

Pi  P2  •  •  •  Pn 

p\  p\  ■■■  Pn 

p r1  pr1  •••  PT 

B.  Parameter  Estimation 

Consider  the  m  x  1  data  vector  {y(t)}  as  given  by  Equation  (1).  In  general  we 
first  filter  the  data  set  {i/(t)}  before  decimating,  to  minimize  effects  of  aliasing  (we 
discuss  filtering  in  detail  in  Section  IV).  Here,  we  consider  an  /th  order  FIR  filter  of 
the  form 


y\t)  =  Ky(t ), 


(5) 


where 


k\  ki- 1  &/_  2  &o  0  0  •  •  •  0 

0  ki  ki-y  ■  •  •  0  •  •  ■  0 

:  ;  ,  (6) 

0  0  ki  fc/_i  •  •  •  ki  k0  0 

0  •  •  •  0  0  ki  ■■■  k2  h  k0 

-■  (m—l)xm 

and  where  the  sequence  {^c)c=o  is  the  FIR  filter  impulse  response. 

From  y'q(t)  we  now  define  a  set  of  decimated  sequences  as 

y?(i)  =  y'qd+uV)  9  =  0,1,...,  m'd  —  1  u  =  0,1,. ..,d-  1,  (7) 

where  m'd  =  L^^J-  The  index  “u”  gives  the  start  sample  in  the  decimation;  thus, 
the  sequences  (for  fixed  t)  represent  a  set  of  interleaved  sequences  decimated 

from  {y'{t)}.  From  Equations  (2)  and  (7),  we  see  that  each  sequence  {y'qn(t)}  is  a 
noisy  exponential  sequence  of  the  form 

y’u(t)  =  ±xnt)W',  +  e?(t),  (8) 

i= 1 

where 

P\  =  (p.y 

x?W  =  x,(()p!“+0K(p,).  (9) 

and  where  £(2)  is  the  FIR  filter  polynomial  given  by 

JC(z)  —  ko  kjz  1  +  k2z  2  +  •  •  •  -f  kiz  (10) 
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Note  that  the  effects  of  the  FIR  filtering  in  the  new  model  are  to  scale  the  ampli¬ 
tude  coefficients  and  to  color  the  noise.  In  general,  we  will  choose  K  to  be  a  bandpass 
filter.  By  careful  choice  of  the  FIR  filter,  we  can  significantly  reduce  the  mode  am¬ 
plitude  coefficients  outside  of  some  band  of  interest;  in  this  case  we  can  assume  the 
number  of  the  “significant”  modes  in  the  filtered  data  is  n'  which  is  less  than  n.  In 
this  case,  we  have  the  following  model: 

jC(‘>  =  (ii) 

1=1 

where 

<*«=  £  *?(t)  (?;)’  + <“(<)•  (12) 

i=n'+ 1 

Note  that  e'“(f)  is  colored  Gaussian  noise  with  nonzero  mean;  the  effect  of  the  nonzero 
mean  is  to  introduce  some  bias  in  the  parameter  estimates,  as  we  will  see  in  the 
following  sections. 

Thus  we  now  have  N  x  d  sequences  with  common  poles  but  different  amplitude 
coefficients.  This  case  is  similar  to  the  original  multi-snapshot  nondecimated  model 
in  Equation  (2).  As  a  result,  the  TLS-Prony  algorithm  can  be  applied  to  the  data 
in  Equation  (7)  to  give  estimates  and  j.  We  thus  have  decimated  multi¬ 

snapshot  backward  linear  prediction  equations  given  by: 

y'°  Y'° 

yn  y/I 

y>  d~l  Y,d~  1 
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where 


>  H [£] 


and  each 


b'  =  b\  b'2  ■■■  b'L 


J/u+dU) 


y'u+d(l) 

J/u  +  2d(1) 


^u  +  2d(l) 

2/u  +  3d(^) 


y'u+LdW 

24+(L+l)d(1) 


y'u  y, 


—  (C  +  l  )d^  )  3/tt+m^  —  Ld(  ^  — (L  —  l)d^) _  —  d(~0 

24(2)  3/'+d(  2)  24+2d(2)  yi+id(2) 

y'+d(2)  2/L  +  2d(2)  2/„+3d(2)  •  •  yu+(X,+l)d(2) 

y’u  +  m^-(L  +  l)d(2)  2/(i+m^-Ld(2)  2/u+m^-(/,-l)d(2)  2/u+m^-d(2) 


y(,(Ar) 

yUd(Ar) 


2 Ud(N) 

y'+2d(Af) 


y'+2d(Ar) 

24+3d(^) 


•••  y'»+Ld(N) 

2/„+(i+i)d(Ar) 


y„+m^_(i,+i)<i(A'’)  y'u+m't-Ld^)  y'u+m'd-(L-i)d(N)  ■■■  yu+m;,-d(AO 


Here  I  is  the  order  of  prediction,  and  b'  is  the  coefficient  vector  of  the  polynomial 
B\z)  given  by 

B\z)  =  1  +  b\z  +  b'2z2  +  •  •  •  +  b'LzL.  (17) 


The  choice  of  L  affects  the  accuracy  of  the  b\  coefficients,  as  we  address  in  later 
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sections. 


The  TLS-Prony  method  considers  the  effect  of  noise  perturbation  of  both  Y'  and 
y',  and  the  TLS  solution  attempts  to  minimize  the  effect  of  these  perturbations  on  the 
prediction  coefficient  vector  b'  (see  [1]  for  details).  This  is  accomplished  by  obtaining 
an  SVD  of  the  matrix  y'  ;  y7  and  truncating  all  but  the  first  n'  singular  values 
to  arrive  at  an  estimate  y>  :  Y'  [lj.  Inserting  y>  :  Y'  in  Equation  (14)  gives 
the  modified  linear  prediction  equation 

Y'b'  ss  -y'  (18) 

from  which  the  linear  prediction  coefficient  vector  estimate  b'  is  found  as 

v  =  - Y'+y (19) 

where  +  denotes  the  Moore-Penrose  pseudoinverse.  Finally,  the  estimates  for  the 
decimated  poles  are  found  by 

p'j  =  zero,-  (£'(*))  ,  j  =  1, 2, . . . ,  L.  (20) 

It  is  not  in  general  possible  to  recover  pj  from  p'  ,  as  the  mapping  p,  — +  p^  is  not 
one-to-one.  However,  the  mapping  can  be  made  one-to-one  by  suitably  restricting 
the  domain  of  pt.  For  example,  if  it  is  known  a  priori  that  Zp,  €  then  Pj 

may  be  uniquely  recovered  from  p'  by 

Pi  =  {p'j)d-  (2!) 
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In  order  to  meet  the  domain  restriction  requirement  without  a  priori  information, 
one  can  choose  a  suitable  FIR  filter  K ,  as  discussed  in  Section  IV  below. 

Once  the  poles  are  found,  the  corresponding  amplitude  coefficients  can  be  esti¬ 
mated  from  the  decimated  pole  estimates  and  decimated  data  using  Equations  (11) 
and  (9)  or  from  the  nondecimated  pole  estimates  and  nondecimated  data  in  precisely 
the  same  way  as  in  the  nondecimated  TLS-Prony  solution  [10].  Using  the  decimated 
pole  estimates  and  decimated  data  is  more  computationally  efficient,  since  there  will 
be  shorter  data  lengths  and  fewer  pole  estimates.  For  this  case,  Equation  (11)  leads 
to  the  following  equation  for  the  amplitude  coefficients, 

A'lKpX  =  Y:°,  (22) 

where 

1  1  •••  1 

Pi  P2  ■■■  Pl 

P\  ?2  ■■■  Pl 

Pi  P2  •••  Pl 

KP  =  diag  (pl'lCifi)  ,...,pl‘)C{pl)'j 

xj(l)  £j(2)  •••  Xi(N) 

2:2(1)  ^2(2)  •••  x2(N) 

xL(l)  xl(‘2)  ■■■  xl{N) 
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The  amplitude  coefficients  can  be  found  from  a  least  squares  solution  to  Equation  (22), 

X  =  £-' (A'lA'LyJ  A'lY«,  (24) 

where  *  denotes  complex  conjugate  transpose.  We  note  that  Equation  (24)  is  not 
used  in  practice  to  solve  Equation  (23),  as  more  numerically  sound  procedures  (such 
as  QR  decomposition  [11])  can  be  used. 

Because  only  n'  singular  values  of  Y'  are  nonzero,  there  are  at  most  n'  pole  esti¬ 
mates  which  can  correspond  to  true  data  modes.  Therefore,  only  the  n'  poles  which 
have  the  largest  energy  are  retained  [10].  We  then  re-estimate  the  amplitude  coef¬ 
ficients  of  these  n'  poles  by  eliminating  all  but  the  n'  “high  energy  pole”  columns 
of  A'l,  then  recomputing  the  least  squares  solution  for  X.  We  note  that  the  second 
amplitude  coefficient  estimation  can  be  done  by  using  the  QR  decomposition  from  the 
first  amplitude  coefficient  estimation.  By  doing  so  we  save  computation  in  the  second 
amplitude  coefficient  estimation.  We  also  note  that  the  noise  e'qu(t)  in  Equation  (11) 
is  not  in  general  white  so  an  (unweighted)  least  solution  to  Equation  (22)  may  not 
lead  to  amplitude  estimates  with  smallest  variance. 

The  above  procedure  can  be  modified  further  to  provide  even  greater  computa¬ 
tional  savings.  The  extra  interleaving  decimated  data  sets  can  be  discarded  for  the 
computation  of  the  poles,  i.e.,  u  in  Equation  (7)  can  take  on  only  the  value  0.  From 
a  Nyquist  theory  point  of  view,  we  note  that  discarding  the  interleaving  data  should 
provide  no  loss  in  performance  if  an  ideal  lowpass  filter  can  be  implemented  (this  is 


confirmed  in  the  examples  presented  in  Section  V  below).  Note  that  this  data  dis¬ 
carding  can  be  easily  incorporated  into  the  matrix  K  by  keeping  only  every  dth  row 
and  that  the  statistical  analysis  below  is  presented  in  a  general  K  framework,  so  that 
the  analysis  applies  equally  well  to  the  cases  where  data  is  and  is  not  discarded. 


III.  Statistical  Analysis 


A  major  contribution  of  this  paper  is  the  derivation  of  the  statistical  proper¬ 
ties  of  the  TLS-Prony  pole  estimates  obtained  using  decimation.  Below  we  derive 
a  general  expression  for  the  first  order  approximation  of  the  probability  distribution 
function  (pdf)  of  the  estimates  of  {p'}  under  the  assumptions  that  there  is  a  filter 
present  as  described  by  Equation  (5).  This  expression  applies  to  different  decima¬ 
tion  values  so  it  can  be  used  to  determine  the  relative  statistical  accuracy  for  various 
choices  in  the  TLS-Prony  algorithm.  The  expression  is  given  in  the  following  theorem. 

Theorem  4-1:  Assume  we  are  given  FIR  filtered  data  {y'(t)}  as  defined  in  Equa¬ 


tions  (7)  and  (11).  Let 


P'  ~  P'i  P'i 


be  the  n'  highest  energy  TLS-Prony  pole  estimates  found  from  Equations  (20)  and 
(24).  Then  the  first  order  approximation  (as  a  —*  0)  of  the  pdf  of  P'  is  given  by 


P'  ~  M  (P'  +  Pl,Z) ,  (26) 
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where  the  estimate  bias  and  covariance  are  given  by 


Pi  =  F'G'S,+Pl , 


£'  =  aF'G'S,+Z £  (F'G'S'+y 


Here  P[  is  defined  in  Equation  (68)  and  £t  is  given  by 


£°  £3  •  •  •  YA 


£}*  £°  ••• 


£rf-!*  •••  £J*  £° 


where  each  £[v  “1  is  a  block  diagonal  matrix  given  by 


=  diag  (B'1<UKV’B'\  B'KuKv'B”,  . . . ,  B'KuKvmB'*) 


N(m'd-L)xN(m'-L) 


The  expressions  for  Ku  and  B'  are  given  by 


1 T 

=  K  (u  +  1  )J  K(u  +  ]+d)T  ■■■  I\  (u  +  1  +  m!  —  d)T  >  (30) 


where  K(a)  is  the  oth  row  of  K,  and 


1  b[  b'2 


b',  0 


0  1  b\  ■■■  b'L_,  b'L 


O'"  o  1  b\  ...  b'L_ ,  b'L  0 


0  0  0  1 


h- 2  hL- 1  VL 


(m'-L)x(m') 
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(32) 


(33) 


(34) 


Proof:  See  Appendix.  □ 

The  above  theorem  gives  the  accuracy  expressions  for  the  case  where  the  data  are 
decimated  by  a  factor  of  d,  and  all  d  interleaved  data  sets  are  used  in  the  estimation 
(see  Equation  (13)).  In  some  cases  one  may  wish  to  use  only  one  of  the  d  interleaved 
data  sets,  in  which  case  Equation  (13)  is  replaced  by 


where  y'°  and  V''0  are  defined  in  Equation  (16).  For  this  case,  the  accuracy  expressions 
for  the  resulting  parameter  estimates  are  given  in  the  following  corollary. 

Corollary  f.l:  Assume  we  are  given  data  as  in  Theorem  f.l,  but  form  TLS-Prony 
estimates  based  on  Equation  (35)  instead  of  Equation  (13).  Then  the  first  order 
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(36) 


approximation  {as  o  -*  0)  of  the  pdf  of  P'  is  given  by 

where 


P'bl  =  F'G'S,0+  P®, 

s;  =  aF'G'S,0+T,°c  (F'G'S,0+y.  (37) 

Here  P'c°  is  the  first  block  of  P[  (the  block  involves  A'0),  and  E°,  F',  and  G'  are  the 
same  as  the  ones  in  Theorem  f.l,  and  S'°  is  the  noise  free  version  of  Y'°. 

Proof:  Follow  by  replacing  5"  by  S'°  and  I<  by  K°  in  the  proof  of  Theorem  f.l.  □ 
Equations  (27)  and  (37)  provide  the  biases  and  covariances  for  decimated  pole 
estimates  given  a  particular  set  of  poles  and  decimation  factor.  If  the  nondecimated 
pole  estimates  {pj}  are  recovered  from  the  decimated  pole  estimates  {p^}  using  Equa¬ 
tion  (21),  the  biases  and  variances  for  the  nondecimated  poles  can  be  derived  in  terms 
of  the  biases  and  variances  for  the  decimated  poles  using  a  first  order  approximation 
of  Equation  (21).  Defining  p-  =  p(  —  p(.  We  have  the  following  derivation. 

Pi+Pi  =  {Pi)d 

=  '  (Pi+PzY 

—  Pi  +  +  (higher  order  terms),  (38) 
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Therefore,  we  have 


Bias(pi)  = 
Var(pt)  = 


E  [pi]  -  pi  = 


Bias  (p[) 

dpt 1 


E  [(pi  -  E  [pi])  (p,  -  £  [pi])*]  = 


Var  (p-) 

d 2  |pi|2(a£_1)' 


(40) 


The  variances  of  the  nondecimated  poles  can  now  be  compared  to  their  respective 
CRBs  (provided  the  estimated  bias  is  negligible).  Recently,  a  CRB  formulation  has 
been  developed  for  multi-snapshot  damped  exponentials  in  noise  [12].  These  CRB 
results  can  be  directly  compared  with  the  variances  of  the  estimated  poles  using  the 


TLS-Prony  method  to  examine  its  performance  in  both  nondecimated  and  decimated 
circumstances.  This  comparison  is  shown  in  Section  V  for  a  number  of  examples. 


IV.  Full  Spectrum  Estimation  Using  Filtering  and  Decimation 

Using  the  decimation  scheme  which  has  been  developed,  any  poles  or  modes  not 
in  the  band  of  interest  need  to  be  filtered  out  so  that  they  are  not  aliased  into  the 
band  of  interest  by  the  decimation  operation.  Even  if  there  are  no  poles  outside  the 
band  of  interest,  a  filter  can  be  still  applied  to  reduce  the  out-of-band  noise  which 
will  be  aliased  into  the  band  of  interest  by  the  decimation  operation.  However,  the 
imperfections  of  a  FIR  filter,  nonideal  stopband  rejection  and  data  length  reduction  by 
transient  response  effects,  will  cause  loss  in  performance,  as  shown  in  Section  V  below. 
In  general,  nonideal  stopband  rejection  increases  bias  in  the  estimation  (because  the 
leakage  of  the  stopband  poles,  and  thus  the  first  term  in  Equation  (12),  will  be  larger), 
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and  data  length  reduction  due  to  transient  effects  of  the  FIR  filter  causes  a  variance 
increase. 

In  this  section,  we  develop  a  procedure  to  obtain  full  spectrum  mode  estimates 
by  use  of  bandwidth  segmentation,  and  then  discuss  how  to  design  the  needed  FIR 
filters  so  that  the  performance  loss  is  minimized.  Finally,  we  present  operation  counts 
for  both  the  decimated  and  nondecimated  estimation  procedures. 

A.  Procedure 

To  examine  the  use  of  filtering  and  decimation,  assume  we  are  interested  in  es¬ 
timating  modes  in  the  complete  spectrum  (that  is,  poles  which  may  lie  anywhere  in 
the  complex  plane),  and  we  wish  to  use  a  decimation  factor  of  d.  We  estimate  poles 
in  each  region  of  Figure  1(a)  using  a  decimation  based  TLS-Prony  procedure.  First, 
we  modulate  the  data  to  center  the  band  of  interest  about  f  —  0  as  follows 

y9(t)  ->  y<,{i)e-:2*Ioq,  (41) 

where  /0  is  the  modulating  frequency  for  one  of  the  subbands  of  interest,  f0  = 
0,  q, . . . ,  We  then  lowpass  filter  the  modulated  data  to  isolate  the  frequency 
band  /  £  ~u'2d}-  Finally,  we  apply  the  decimated  TLS-Prony  algorithm  of  Sec¬ 
tion  II.  The  resulting  pole  estimates,  p',  as  given  by  Equation  (21)  lie  in  /  £  — ~ 
as  shown  in  Figure  1(b).  The  corresponding  pole  estimates  in  nondecimated  frequency 
space  are  given  by  Equation  (21)  with  modulation  as  follows 

Pi  =  (p'j)d  ej2*fo-  (42) 
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(a)  (b) 


(c) 

Figure  1:  Spectral  effect  due  to  decimating  by  six. 

A  problem  which  can  result  from  the  above  procedure  is  that  poles  near  the 
endpoints  of  the  subband  region  may  be  incorrectly  estimated  near  the  opposite 
endpoint.  This  results  from  the  discontinuity  of  the  mapping  in  Equation  (42)  for 
Zp(  «  ±7r.  It  can  be  seen  from  Figure  1(b)  that  small  errors  in  estimates  near 
A  and  B  result  in  large  differences  in  Figure  1(a).  To  avoid  this  problem,  we  use 
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c  >  d  overlapped  estimations,  each  of  size  ^  (note  that  this  changes  the  modulating 
frequencies  to  be  /0  =  0,  |, . . . ,  ^).  This  is  shown  in  Figure  1(c)  for  c  =  2d.  For 

each  region,  we  retain  only  those  pole  estimates  which  are  in  the  half  of  the  overlap 
region  which  is  closer  to  the  center  of  the  band  of  interest.  For  example,  in  Figure  1(c) 
for  the  subband  [A,B]  we  retain  poles  only  in  the  region  [C,D].  This  corresponds  to 
retaining  pole  estimates  whose  angles  satisfy 


in  the  decimated  frequency  space.  In  this  way,  we  reduce  the  effects  of  the  disconti¬ 
nuity  of  the  mapping  in  Equation  (21)  for  Ip]  near  n.  The  overlap  method  also  helps 
to  provide  immunity  to  effects  of  a  nonideal  lowpass  filter,  as  is  discussed  in  the  next 
subsection. 


B.  Filter  Design  and  Performance  Loss 

As  we  discussed  before,  the  use  of  finite  length  FIR  filters  results  in  performance 
loss.  We  need  to  design  a  filter  such  that  this  performance  loss  is  minimized.  Because 
the  stopband  of  the  filter  is  set  by  the  decimation  factor,  d,  to  be  there 

remain  three  free  design  parameters:  the  order  of  the  FIR  filter  (/),  the  number  of 
overlapped  estimations  (c),  and  the  flatness  of  the  passband.  Note  that  the  passband 
is  defined  as  —  T  ^  the  region  in  which  we  actually  retain  the  pole  estimates. 

Each  of  the  parameters  has  its  own  effect  in  the  estimation  performance.  Larger 
filter  lengths  result  in  a  variance  increase  because  of  the  data  length  reduction  associ- 
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ated  with  discarding  the  transient  response  of  the  filter  output.  Shorter  filter  lengths 
result  in  increased  bias  of  the  pole  estimates  due  to  aliasing  of  imperfectly  attenu¬ 
ated  modes  in  the  stopband  region.  The  number  of  overlapped  estimations  primarily 
influences  the  computational  aspect.  A  larger  c  results  in  more  computations;  never¬ 
theless,  a  larger  c  allows  for  a  larger  filter  transition  band,  thus  allowing  us  to  design 
for  a  better  filter  stopband.  Note  that  we  can  tolerate  some  non-flatness  since  we  can 
compensate  for  this  effect  in  the  amplitude  coefficient  estimation.  However,  a  non-flat 
passband  increases  the  variance  of  the  poles  near  the  minima  of  the  passband  since 
the  filter  reduces  the  power  of  these  poles  more  than  other  poles  in  the  passband. 

After  filtering  and  decimating  the  data,  the  stopband  signal  poles  are  aliases  into 
the  band  of  interest.  If  the  energy  of  the  filtered  stopband  signal  modes  is  small 
relative  to  the  in-band  Gaussian  noise,  the  estimation  bias  caused  by  the  stopband 
signal  poles  will  be  negligible  in  comparison  to  the  estimation  variance.  This  implies 
that  for  high  SNR  signals  we  need  more  stopband  rejection  in  order  to  avoid  the  bias 
problem.  One  should  choose  a  filter  length  sufficient  to  attenuate  out-of-band  modes 
to  below  the  noise  floor. 

A  procedure  to  design  the  FIR  filter  is  presented  here  to  minimize  variance  and 
bias.  First,  based  on  SNR  one  can  determine  the  needed  stopband  rejection  to  ensure 
the  leakage  of  the  out-of-band  signal  poles  is  small  relative  to  the  in-band  noise. 
Then,  one  can  choose  the  filter  order  l  and  the  number  of  the  overlapped  estimations 
c  to  give  the  desired  stopband  rejection  with  a  small  transition  band.  Note  that  we 
want  l  and  c  to  be  as  small  as  possible.  For  instance,  in  the  example  given  below,  we 
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estimated  according  to  the  SNR  that  a  20dB  stopband  rejection  is  enough  to  avoid 
the  bias  problem.  After  few  trials,  we  found  that  using  c  =  2d  =  12  and  1  =  20, 
an  equiripple  FIR  filter  gives  the  approximated  stopband  rejection.  Compromise  was 
made  to  sacrifice  the  flatness  of  the  passband  a  little  since  we  can  compensate  its 
effect  during  the  amplitude  coefficient  estimation. 

C.  Operation  Counts 

The  main  goal  of  applying  the  decimation  procedure  is  a  reduction  in  computa¬ 
tions.  In  this  section  we  compare  the  computations  of  a  full  TLS-Prony  estimation 
applied  to  a  non  decimated  sequence  to  a  set  of  those  from  a  set  of  TLS-Prony  esti¬ 
mates  obtained  from  a  set  of  bandpass  filtered  and  decimated  sequences.  We  include 
the  computations  associated  with  overlapping  as  well  as  those  associated  with  the 
filtering  operation  for  this  comparison. 

In  the  operation  counts  which  follow,  we  assume  that  only  one  of  the  interleaved 
sets  of  decimated  data  is  used  in  the  decimation-based  TLS-Prony  algorithm  (as  in 
Equation  (35)).  In  addition,  we  assume  Z,  =  ^  is  used  for  the  prediction  order, 
because  this  prediction  order  gives  near  optimal  accuracy  (see  [10,  13]  and  Section  V 
below).  We  compute  operation  accounts  for  the  SVD  operation,  the  QR  decompo¬ 
sition,  the  polynomial  root  finding  operation,  and  the  filtering  operation;  these  were 
found  to  account  for  over  90%  of  the  total  operations  in  our  computer  simulations. 
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i.  SVD  Operation  Counts 


For  a  real  matrix  which  has  dimension  R  x  C,  the  approximate  floating  point 
operation  (flop)  count  associated  with  the  “economical”  SVD  computation  (in  which 
only  the  first  C  left  singular  vectors  are  computed)  is  given  by  fc^^  «  l ARC2  + 
8C3  [11].  For  a  complex  matrix  the  count  is  about  a  factor  of  2  larger.  In  our  case, 
the  matrix  y>  y>  has  dimension  R  x  C  —  x  (^57^  +  l)-  Therefore, 

for  the  nondecimated  case  (d  =  1  and  /  =  0)  we  thus  get  an  SVD  flop  count  of 

t4V„Ddec^8(|mW)(9  +  l)2+,6(=  +  l)3.  (44) 


and  for  the  decimated  case  with  c  overlap  regions  we  obtain 


1  +  16 


m  —  l 


2.  QR  Decomposition  Operation  Counts 

For  a  real  matrix,  the  approximate  flop  count  associated  with  the  QR  decompo¬ 
sition  is  given  by  fcQ^  ~  IRC2  —  |C3  [11].  For  complex  matrices,  the  count  is  to  be 
a  factor  of  4  larger.  In  our  case,  R  x  C  =  x  (^3^)-  We  thus  obtain  for  the 

nondecimated  and  decimated  cases,  respectively 

j-cQR  ^  8m3  _  8m3 
1  iiondec  ~  9  gj  ’ 

and 

f  QR  ^  -  ( 8(??l  -  03  S{m-l)3\ 

,cdec~c^  9ds  ~  81  c?3  )■ 
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3.  Polynomial  Root  Finding  Operation  Counts 


For  a  real  polynomial  of  order  k,  the  approximate  flop  count  associated  with  the 
root  finding  operation  is  given  by  fcroot  m  f  /c3  [11].  For  complex  data  considered 
the  count  is  to  be  a  factor  of  about  10  larger;  in  our  case  k  =  So  we  have 

r  root  200  3 

cnondec  K  gj  m  >  (48) 


and 


r  root 
lcdec 


(49) 


4.  Filtering  Operation  Counts 

For  the  decimated  case,  we  must  also  include  the  operation  count  associated  with 
the  filtering  operation.  From  Equations  (5)  and  (6),  the  approximate  flop  count 
associated  with  the  matrix  multiplication  induced  by  this  filter  is  thus  given  by 

fcd^Cr«c(!^(/+l))  x6-  (50) 

The  factor  of  6  arises  because  there  are  6  flops  (4  multiplies  and  2  adds)  per  complex 
multiplication. 


V.  Examples  and  Simulation  Studies 

Examples  using  the  statistical  analysis  results  are  presented  here  which  demon¬ 
strate  the  advantages  of  using  decimation.  Simulations  are  also  presented  for  full 
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spectrum  data  sets  to  demonstrate  the  estimation  ability  of  the  modified  TLS-Prony 
method  developed  in  Section  IV.  Again,  note  that  we  are  considering  the  case  where 
the  extra  interleaved  data  sets  are  discarded. 

A.  Single  Undamped  Mode 

In  this  example  we  assume  one  snapshot  of  data  of  length  m  =  140.  We  assume 
a  single  exponential  located  on  the  unit  circle  and  an  SNR  of  5dB.  We  compare  esti¬ 
mates  of  this  exponential  using  no  decimation  and  using  decimation  by  d  =  6.  Using 
the  filter  design  procedure  outlined  above,  we  obtain  an  FIR  filter  whose  frequency 
response  is  shown  in  Figure  2.  Here  we  used  c  =  2d  overlapped  estimations  and  an 
equiripple  FIR  filter  of  order  /  =  20.  Note  that  we  could  decrease  the  order  of  the 
FIR  filter  to  achieve  the  same  stopband  rejection,  but  in  doing  so  we  will  obtain  a 
less  flat  passband,  which  results  in  increased  variance.  Note  also  that  in  this  case  we 
do  not  have  out-of-band  signal  poles,  so  there  will  be  no  bias  in  the  pole  estimates 
(to  a  first  order). 

Figure  3  shows  the  theoretical  variance  of  the  estimated  pole  versus  prediction 
order  for  various  decimation  factors  as  compared  to  the  CRB.  From  this  figure  we 
can  see  that  the  minimum  variance  occurs  at  a  lower  prediction  order  for  higher  d. 
The  minimum  occurs  for  a  prediction  order  equal  to  about  one  third  of  the  decimated 
data  length  (be.,  L  =  —*),  which  is  consistent  with  results  for  nondecimated  data  [10, 
13].  This  shows  that  the  best  performance  for  the  decimated  cases  occurs  at  lower 
prediction  orders  than  for  the  nondecimated  case,  thus  reducing  the  computational 
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Magnitude  response  of  the  20th  order  FIR  filler 


Figure  2:  Frequency  response  of  the  20 th  order  equiripple  FIR  filter. 

load. 

We  note,  however,  that  since  the  data  has  to  be  filtered  prior  to  decimation,  the 
curves  for  the  decimated  cases  peak  at  about  1.5dB  lower  than  for  the  nondecimated 
case  for  various  decimation  factors.  The  performance  loss  is  due  to  the  fact  that 
the  transient  response  portion  of  the  filter  output  (20  points  for  this  case)  needs  to 
be  discarded.  Note  that  for  a  fixed  length  FIR  filter  this  performance  degradation 
becomes  smaller  as  the  data  length  is  increased,  since  the  percent  difference  between 
the  original  and  filtered  data  lengths  decreases. 
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Figure  3:  1 U  log10  (CRB/Var(p))  for  various  d  and  prediction  orders  for  a  single  un¬ 
damped  mode. 

B.  Two  Undamped  Modes 

In  a  second  example  we  make  the  same  assumptions  above,  except  that  there  are 
two  equal  energy  exponentials  located  on  the  unit  circle  one  Fourier  bin  apart  (be., 
A  / '  =  —  =  T7- ).  The  total  SNR  is  assumed  to  be  8dB  in  this  case  in  order  to  maintain 
5dB  SNR/pole.  Figure  4  shows  the  theoretical  variance  of  the  estimates  for  one  of  the 
poles  (the  variance  of  the  other  pole  is  similar).  We  can  see  that  the  characteristics 
are  much  the  same  as  in  the  one  pole  case,  the  difference  being  higher  variances  due 
to  the  presence  of  each  pole’s  neighbor. 
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Figure  4:  101og10  (CRB/Var(pi))  for  various  d  and  prediction  orders  for  two  un¬ 
damped  modes. 

C.  Monte-Carlo  Simulation  of  an  Undamped  Ten  Mode  Case 

We  now  present  a  set  of  simulations  for  a  general  ten  undamped  mode  case.  In 
these  simulations  we  have  N  =  1  snapshot,  and  n-  10  poles  present  in  the  data.  The 
amplitude  coefficients  all  have  unit  magnitude;  the  phases  of  the  amplitude  coefficients 
are  chosen  randomly.  We  consider  two  data  lengths,  m  =  140  and  m  =  560  data 
points.  Figure  5  shows  the  locations  of  the  ten  poles;  each  is  indicated  by  an  “x” . 
Five-hundred  independent  Monte-Carlo  simulations  are  performed  by  adding  noise  to 
the  data  such  that  the  total  SNR  is  20dB  (lOdB/pole).  Estimates  for  the  poles  are 
obtained  using  the  TLS-Prony  algorithm  without  decimation  (i.e.,  d  =  1)  and  with 
decimation  using  a  decimation  factor  of  d  =  6. 
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Figure  5:  True  pole  locations  for  an  undamped  ten  mode  case. 

For  the  decimation  results,  the  FIR  filter  is  the  same  as  the  one  used  in  the  previous 
examples  and  the  frequency  response  was  shown  in  Figure  2  (thus,  c  —  12  and  /  =  20). 
The  prediction  order  is  L  —  |_^^J ;  the  numbers  of  singular  values  retained  in  the 
simulations  are  10  for  d  =  1  and  the  number  of  poles  in  each  subestimation  section 
for  d  =  6  ({2,  3, 2, 0, 1, 1, 2, 2, 0,1, 3, 3}  for  this  case).  The  prediction  orders  used 
correspond  to  one  third  of  the  effective  data  lengths  in  the  two  cases  as  was  suggested 
by  Examples  1  and  2.  Prior  to  the  calculation  of  the  amplitude  coefficients,  poles 
with  magnitude  larger  than  1.15  and  smaller  than  ^  are  eliminated  to  avoid  poor 
conditioning  in  the  least  squares  solution  of  the  amplitude  coefficients. 
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Table  1:  Theoretical  and  simulation  variances,  and  MSEs  for  the  undamped  ten  poles, 
in  =  140  data  point  case  (all  values  are  in  dB). 


pole 

d  = 

1 

d  =  6 

number 

CRB 

Theory 

Sim. 

Sim.  MSE 

Theory 

Sim. 

Theory  MSE 

Sim.  MSE 

1 

-63.4 

-61.8 

-61.5 

-61.5 

-60.5 

-60.6 

-60.3 

-60.4 

2 

-63.0 

-60.7 

-60.9 

-60.8 

-59.7 

-58.9 

-59.5 

-58.7 

3 

-63.0 

-60.7 

-61.0 

-60.9 

-60.4 

-58.7 

-59.9 

-58.0 

4 

-63.6 

-61.4 

-61.2 

-61.2 

-60.9 

-61.0 

-60.9 

-60.9 

5 

-63.4 

-62.3 

-62.3 

-62.2 

-60.6 

-60.4 

-60.1 

-59.9 

6 

-63.4 

-62.3 

-62.2 

-62.2 

-60.7 

-60.5 

-60.4 

-60.3 

7 

-63.5 

-61.4 

-61.1 

-61.1 

-61.0 

-60.8 

-60.9 

-60.7 

8 

-55.1 

-51.4 

-51.0 

-50.9 

-47.3 

-48.4 

-47.0 

-48.1 

9 

-55.1 

-51.4 

-51.3 

-51.1 

-46.8 

-47.9 

-46.5 

-47.4 

10 

-63.4 

-61.8 

-61.6 

-61.6 

-59.6 

-59.6 

-59.5 

-59.6 

1.  Performance  Comparison 

Tables  1-3  summarize  the  performance  of  the  various  methods  for  this  example. 

We  first  consider  Table  1  which  shows  the  m  =  140  data  point  case.  From  Table  1 
we  see  that  the  Monte-Carlo  variances  for  the  d  =  1  case  were  1.2dB  to  4.1dB  away 
from  theii  CRBs.  Note  that  the  Monte-Carlo  variances  are  within  0.5dB  of  those 
predicted  by  the  theory  which  substantiates  the  theory  (the  differences  are  due  to 
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the  fact  that  theory  is  only  a  first  order  analysis).  Note  also  that  the  estimates  have 
negligible  bias  which  is  shown  by  the  fact  that  the  pole  variances  are  very  close  to 
their  MSEs. 

The  Monte-Carlo  variances  for  the  d  =  6  case  were  2.6dB  to  7.2dB  away  from 
their  CRBs.  This  represents  an  average  1.7dB  loss  for  the  entire  d  =  6  system 
versus  the  nondecimated  system.  The  main  cause  of  the  performance  loss  is  the  20 
data  point  loss  due  to  the  transient  response  of  the  FIR  filter  outputs.  Note  that 
the  Monte-Carlo  variances  are  within  1.7dB  of  statistical  theory.  By  comparing  the 
simulation  variances  and  MSEs,  we  see  that  some  of  pole  estimates  are  slightly  biased. 
When  decimation  is  used,  Theorem  4-1  gives  an  analytical  expression  for  the  this  bias; 
theoretical  biases  are  compared  to  biases  obtained  from  Monte-Carlo  simulations  with 
good  agreement  in  most  cases. 

11  the  data  length  is  increased  to  m  =  560  points,  the  theoretical  and  simulation 
results  show  even  better  agreement,  as  is  shown  in  Table  3.  In  this  case,  the  overall  loss 
using  decimation  is  less  than  0.4dB  compared  to  the  nondecimated  case.  In  addition, 
the  simulation  variances  are  within  0.4dB  of  the  theoretically  derived  variances,  and 
the  bias  of  the  pole  estimates  is  significantly  reduced. 

2.  Operation  Count  Comparison 

The  d  =  1  and  d  =  6  estimation  procedures  are  now  compared  on  the  basis  of  their 
computational  costs.  We  only  compare  the  results  for  m  =  140  data  points.  Using 
MATLAB,  the  “economical”  version  of  the  SVD  operation,  the  left  division  opera- 
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Table  2:  Theoretical  and  simulation  biases  for  the  poles  (m  =  140,  unit=  10-3) 


pole 

number 

Theory 

Sim. 

1 

0.078-0. 179i 

0.088-0. 196i 

2 

0.166-0. 033i 

0. 255-0. 037i 

3 

0. 330+0. 072i 

0.449+0. 188i 

4 

0. 042+0. 085i 

0. 005+0. lOOi 

5 

-0. 082-0. 314i 

-0. 074-0. 324i 

6 

0. 086+0. 205i 

0. 018+0. 184i 

7 

0. 076-0. 077i 

0.067-0. 1 36i 

8 

1 .251-0. 1 71i 

0. 847-0. 693i 

9 

-1.121-0.470 

-1.1 34-0. 625i 

10 

-0.095+0. 142i 

-0. 047+0. 021i 

tjon  (using  QR  decomposition  to  solve  least  squares  problems),  and  the  root  finding 
operations  required  an  average  of  16.5Mflops  for  each  of  the  Monte-Carlo  simulations 
for  the  d  —  1  case.  Each  of  the  twelve  d  =  6  SVDs,  QR  decompositions,  polynomial 
root  findings,  and  filtering  operations  required  an  average  of  53.7Kflops,  resulting  in 
a  total  of  644.4Kflops  for  each  Monte-Carlo  simulation.  The  computational  savings 
lot  the  SV  Ds,  QR  decompositions,  polynomial  root  findings,  and  filtering  operations 
in  this  example  using  decimation  was  thus  a  factor  of  about  25.6.  This  compares  with 
a  savings  factor  of  24.0  which  is  predicted  by  Equations  (44)-(50)  for  this  scenario. 
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Table  3:  Theoretical  and  simulation  variances,  and  MSEs  for  the  poles  for  the  560 
data  point  case  (all  values  are  in  dB). 


pole 

d=  1 

d  =  6 

number 

CRB 

Theory 

Sim. 

Sim.  MSE 

Theory 

Sim. 

Theory  MSE 

Sim.  MSE 

1 

-81.6 

-79.9 

-79.6 

-79.6 

-79.7 

-79.8 

-79.7 

-79.8 

2 

-81.6 

-79.6 

-79.7 

-79.7 

-79.4 

-79.4 

-79.4 

-79.3 

3 

m 

-79.7 

-79.3 

-79.4 

-79.3 

-79.4 

4 

... 

-79.6 

-79.7 

-79.7 

-78.9 

-79.0 

5 

-81.6 

-79.8 

-79.8 

-79.8 

-79.5 

-79.2 

-79.4 

-79.1 

6 

-81.6 

-79.8 

-79.7 

-79.7 

-79.4 

-79.4 

-79.1 

-78.9 

7 

-81.7 

-79.9 

-79.9 

-79.8 

-79.7 

-79.4 

-79.6 

-79.3 

8 

-81.6 

-79.9 

-80.0 

-80.0 

-79.6 

-79.4 

-79.0 

-79.0 

9 

-81.6 

-80.0 

-79.5 

-79.4 

-79.5 

-79.4 

10 

-81.7 

-80.0 

-79.3 

-79.7 

-79.3 

-79.7 

Note  that  these  numbers  cannot  be  computed  directly  in  part  due  to  differences  be¬ 
tween  MATLAB’s  SVD,  QR  decomposition,  and  root  finding  algorithms  and  the  ones 
used  for  operation  counts  in  Equations  (44)-(50).  The  average  total  flop  counts  (in¬ 
cluding  all  operations)  for  the  nondecimated  and  decimated  Monte-Carlo  simulations 
were  17.4Mflops  and  698Kflops,  respectively,  to  give  a  savings  factor  of  24.9.  Note 
that  the  four  computational  cost  components  we  have  detailed  make  up  about  92% 
of  the  total  computations.  With  higher  decimation  factors  the  savings  are  even  more 
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substantial. 

VI.  Conclusions 

In  this  paper  we  have  developed  a  TLS-Prony  estimation  algorithm  which  incor¬ 
porates  data  decimation.  We  also  have  developed  a  statistical  analysis  for  estimated 
poles  of  this  algorithm.  We  have  shown  through  examples  using  this  analysis  that 
decimation  provides  a  minimum  variance  for  estimated  poles  that  occurs  at  a  predic¬ 
tion  order  which  is  smaller  than  the  optimal  prediction  order  for  nondecimated  data 
by  a  factor  of  d,  thus  allowing  for  computational  savings.  We  have  shown  that  this 
benefit  is  obtained  at  the  expense  of  pole  variance  performance  due  to  the  filtering 
which  is  required;  this  expense  becomes  smaller  for  longer  data  lengths.  We  have 
also  shown  how  the  modified  TLS-Prony  method  can  be  used  on  full  spectrum  data 
one  band  at  a  time  to  realize  the  computational  savings  in  a  more  general  signal 
framework. 

With  this  decimation  procedure,  we  are  now  able  to  make  a  well-quantified  tradeoff 
of  accuracy  for  computation  when  using  the  TLS-Prony  estimation  procedure. 

Appendix:  Proof  of  Theorem 

From  Equation  (14)  we  can  make  the  following  substitutions 

(S' +  S')  (y  +  i')  =  -  (s' +  s') ,  (51) 
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where  s'  is  the  noise  free  version  of  y',  Y'  =  S'  - 1-  S' ,  y'  =  s'  +  s',  and  ^  =  b'  +  b'. 
We  can  see  that  the  •  terms  are  small  perturbations  for  the  high  SNR  case,  which  is 
assumed.  Multiplying  out  Equation  (51)  and  retaining  only  the  first  order  •  terms 
gives1 

S'b'  +  S'b'  +  S'b'  =  -s'  -  's'.  (52) 

Now  note  that  S'b'  =  —s'  since  they  are  the  noiseless  terms.  Equation  (52)  thus 
becomes 

S'b'  =  -  (s'  +  S'b’)  .  (53) 

Multiplying  both  sides  through  by  S'S,+  and  noting  that  S'S’+S'  =  S',  we  obtain 

S'b'  =  -S'S'+  (s'  +  S'b')  .  (54) 

Let  Y'  —  S'  +  W'  and  y'  =  s'  +  w',  where  W'  and.u;'  are  the  appropriate  noise  matrices 
(i.e.,  of  the  form  given  in  Equation  (14)  and  composed  of  the  noise  sequences  in  Equa¬ 
tion  (12)).  Thus  we  can  see  that  S'  and  s'  are  perturbations  caused  by  W' ,  w' ,  and 
the  SVD  truncation.  By  using  perturbation  analysis  [14]  on  the  matrices  s'  S'  i 
y>  y>  ,  and  y>  Y'  >  it  can  be  shown  that  to  first  order  approximation 

S'*  [  y  S' }  =  S'*  [  w  ]  ■  (55) 

Thus,  Equation  (54)  can  be  written  as 

S'b'  =  -S'S'+e,  (56) 

’Note  that  the  approximation  is  valid  since  the  matrices  y<  •  y>  and  s'  :  S'  have 
the  same  rank. 
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where  t  =  w'  -f  W'b1. 

Observing  the  data  model  and  the  formulation  of  the  S'  matrix,  we  can  write  S,u 


*i(l  )P? 
*i(l  )Pi+d 


X2(l)P2 

xi(l)p2+C 


*n'(l  )Pn 
*n'(l  )PUJJ 


,,,  u+m-(Z.  + l)d  /i\  u+m-(Z,+  l)d 

Zl(l)Pl  x2{l)p2 


*i(2  K 
*i(2)p?+< 


*2(2)P2 

*2(2)pS+< 


S'u  = 


*l(2)P! 


ti+m  —  (L+l)d 


a;2(2)p 


u+m  —  (.L+l)<i 


*n'(2)pj{ 

*n'(2)p“+d 


*n'(2)P^m_(i+1)Ci 


•'•i(A'K+C 


*2(A0p“ 

xo(A’)p“+< 


®n'(A^)Pn 


Zn'(A')p“< 


7W,u+'i 


/  a j  \  ti+m  — (Z»+l)«i  /  »r\  u+m  —  (L+l)d  y  *r\  u+m  —  (L  +  l)<i 

zi(A  )Pi  x2(N)p2  ■■■  xn‘{N)pn, 


S'u  =  //“<7. 


T 

Letting  H  =  /y°r  //1?  ...  Hd~yl  ,  we  get 


5'  =  //O'. 


Equation  (56)  thus  becomes 


HG'b 1  =  -HG'S'+c. 


Now  note  by  definition  that  the  true  and  estimated  Lth  order  characteristic  poly¬ 
nomials  are  B(z)  =  1  +  b[z  +  b'2z 2  -| - f  b'LzL  and  B'{z )  =  1  -f  b^z  +  b'2z 2  -f - f  b'LzL, 

respectively.  Hence  B  (pj)  =  0  and  B'  (p()  =  0. 

We  can  use  a  first  order  Taylor  expansion  to  find  an  expression  for  the  error  in 
the  estimated  pole  locations.  For  each  p(  we  obtain 


0  =  B'(p[) 

=  B'  ( )  +  4z^'(z)\z=p[  (p'i  ~  Pi )  +  (higher  order  terms) 

=  B'  (p')  -  B'  (p')  +  ~B'{z)\z=p[  (p'  -  p')  +  ((higher  order  terms) 

~  i  +  b'lVx  +  —  (i  +  K p'i  +  K p?  +  •  •  •  +  Kp'i') 

+  (b\  +  2 v2p\  +  •  •  •  +  LVjjp^)  (pi  -  p') 

f;  -  *>; }  \  i 


Pi  P?  PiL 


b'2  -  b2 


+  K  b'2  •••  61 


(pi  -  Pi) 


b'L-b'L 


P'i  P?  ■ 


j'L  (b>  -  6')  +  rfi  -  p')  , 


AL-*) 


or.  to  first  order. 


(pi  -  Pi) 


7  Pt  P? 


Thus,  for  all  of  the  n'  true  poles  we  obtain 


P'-P'  =  -F'G'b'. 


b 


Since  H  is  full  rank  (this  can  be  seen  by  noting  that  each  block  of  rows  is  simply 
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a  Vandermonde  matrix  derived  from  distinct  poles  times  a  diagonal  matrix  of  the 
nonzero  amplitude  coefficients),  we  can  multiply  Equation  (60)  by  ( H*H)~l  H*  to 
get 

G'b'  =  —G'S,+e,  (64) 

and,  by  substituting  Equation  (64)  into  Equation  (63),  we  obtain 

P'  -  P'  =  F'G'S'+t.  (65) 

We  now  note  that,  w>  if'  |f|,|  can  be  written  as 

£V°(1)  B'K°{sT{l)  +  e{l)) 

B'e,0{  2)  B'K°(sr(2 )  +  c(2)) 

B'e,0(N)  B'K°(sr(N)  +  e(7V)) 

B'en{  1)  £#A'V(l)  +  e(l)) 

B'enC2)  B'K\sr(2)  +  e(2)) 

w1  w  =  ;  :  ,  (66) 

B'en(N)  B'K1(sT(N)  4-  e(N)) 

B'e'd~'(l)  B'I\d~i(sT(l)  4-  e(l)) 

B'c'd~x{2)  B'Kd-\sr{2)  +  c(2)) 

B'e,d~1(N)  B'Kd~i(sr(N)  +  e{N)) 
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where 


T 

*'«=  n.„.+i w‘  •••  E”„.+1^(<)fer‘' j  •  <67) 

and  where  B'  and  “  are  given  by  Equations  (30)  and  (31).  Recall  (e(<)}  are  zero 
mean  Gaussian  and  thus  t  is  multivariate  Gaussian  with  mean 


E  [e]  = 


where  is  in  the  form  of  Equation  (66)  without  the  e(f)’s,  and  covariance  matrix 

Cov(£)  =  E  ([»'  ][*])([«/  W"][i]) 

=  <rS„  (69 


where  Et  is  given  by  Equation  (28). 

Equations  (65),  (68)  and  (69)  imply  that  the  mean  and  covariance  matrix  of  P' 
are  given  by  Equation  (27).  1=1 
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tion  Results  for  Scattering  Centers,”  in  Proceeding  of  the  International  Radar  92 
Conference,  October  12-13,  1992,  Brighton,  UK,  pp.  518-521. 
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'The  Analytic  Saeocee  Corporation,  USA 
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INTRODUCTION 

This  paper  is  concerned  with  detection  and  estimation  of 
the  scattering  centers  of  a  target  from  coherent,  stepped 
frequency  measurements.  In  particular,  we  are  interested 
in  the  following  questions:  1)  how  closely  spaced  can 
scattering  centers  be  before  it  is  impossible  to  resolve 
them,  and  2)  what  is  the  relationship  between  tbe  detec¬ 
tion  probability  of  a  scattering  center  and  the  false  alarm 
probability  as  a  function  of  scattering  center  SNR. 

lb  address  these  questions,  we  hypothesise  a  paramet¬ 
ric  model  of  target  scattering.  This  model  assumes  tbe 
frequency -domain  scattering  to  be  a  sum  of  exponential 
terms.  If  tbe  exponential  terms  are  undamped,  then  the 
model  specialises  to  a  point -scat  terer  assumption.  If  the 
exponential  terms  are  not  undamped,  the  model  incorpo¬ 
rates  frequency-dependent  radar  cross  section  of  scatter¬ 
ing  centers  (see  (1)).  We  consider  a  particular  class  of  al¬ 
gorithms  for  estimating  the  parameters  in  this  exponen¬ 
tial  model,  the  so- called  total-least  squares  (TLS)  Prony 
algorithm  (2).  The  TLS-Prony  technique  has  gained 
popularity  as  a  parameter  estimation  algorithm  for  the 
exponential  model  because  it  provides  accurate  parame¬ 
ter  estimates  at  moderate  computational  cost  (3);  it  has 
also  been  successfully  applied  to  both  one-dimensional 
and  two-dimensiaonal  radar  scattering  data  (4,5). 

Under  the  exponential  model  assumption  of  scattering, 
the  resolution  and  detection  bounds  can  be  reformu¬ 
lated  in  terms  of  parameter  estimation  accuracy  for 
the  exponential  model.  We  present  the  asymptotic  (as 
SNR  — •  oo )  probability  density  function  (pdf)  for  the 
exponential  model  parameter  estimates  using  the  TLS- 
Prony  algorithm.  We  then  use  this  asymptotic  pdf  to 
derive  scattering  center  resolution  and  detection  bounds 
for  the  TLS-Prony  algorithm,  and  compare  these  results 
to  Cramer-Rao  bound  (CRB)  results.  Monte- Carlo  sim¬ 
ulations  are  also  presented  to  compare  with  the  theory. 
The  resolution  bounds  are  obtained  from  tbe  standard 
deviation  bounds  of  tbe  pole  angles  in  the  exponential 
models.  The  detection  bounds  are  obtained  by  consid¬ 
ering  tbe  probability  that  tbe  energy  of  an  estimated 
mode  exceeds  a  pre-defined  threshold.  In  each  case,  tbe 
probabilities  ore  obtained  by  considering  *  high-SNR  ap¬ 
proximation  of  tbe  statistic*!  probabilities  of  exponential 
mode]  parameter  estimates. 

One  of  the  advantages  of  s  model-based  scattering  cen¬ 
ter  estimation  procedure  is  tbe  capability  of  resolving 
scattering  centers  more  accurately  than  is  possible  using 
fast  Fourier  Transform  (FFTJ  techniques.  We  shov  that 
for  sufficiently  high  SNR,  both  the  CRB  resolution  and 
(he  TLS-Prony  resolution  is  better  than  can  be  obtained 
using  the  FFT. 


•THIS  RESEARCH  WAS  SUPPORTED  IN  PART 
BY  THE  AIR  FORCE  OFFICE  OF  SCIENTIFIC  RE¬ 
SEARCH,  THE  AVIONICS  DIVISION,  WRIGHT  LAB¬ 
ORATORIES,  AND  THE  SURVEILLANCE  DIVISION, 
ROME  LABORATORIES 


DATA  MODEL  AND  TLS-PRONY 
ESTIMATION  PROCEDURE 

Data  Model 

Assume  the  data  vector  p  of  length  m  is  modeled  as  a 
noisy  exponential  sequence 


».  =  £*iPf  +  e,  *  =  0,1 . m—  1.  (1) 

•ml 

There  are  n  distinct  exponential  modes  in  the  data. 
Here,  it  is  assumed  that  {et}  is  a  aero  mean  complex 
white  Gaussian  noise  sequence  with  variance  o .  Equa¬ 
tion  1  may  be  compactly  written  as 

»  =  A*  +  e,  (2) 

where  e  (m  x  1)  and  x  (n  x  1)  are  vectors,  and 


A  ~ 


1 

R 


1 

Pi 


1 

P- 


Pi  Pa 


(3) 


TLS-Prony  Estimation  Procedure 

In  this  subsection  we  give  on  overview  of  the  TLS-Prony 
technique  (2)  which  ii  need  to  estimate  the  parameters 
of  th«  exponential  model  presented  in  Equation  2. 


Pint,  Lxh  order  backward  linear  prediction  equations  are 
formed: 


l»  Y] 

!] 

fc  0, 

(4) 

where 

h  =  !  h. 

... 

h  f 

(5) 

and  where 

l»  r  1  = 

'  IS 

IS 

rfi.  ' 

-  (6) 

■  Pm- 1  J 


In  genera],  L  >  n;  however,  choosing  L  >  n  results  in 
more  accurate  parameter  estimates  (6). 

The  solution  of  Equation  4  involves  obtaining  an  S VD  of 
tbe  matrix  [  y  Y  )  and  truncating  all  but  tbe  first  n 
singular  values  to  arrive  at  an  estimate  [  y  Y  ].  Tbe 
linear  prediction  coefficient  vector  estimate  6  is  found  as 
5  =  —  y4y,  where  4  denotes  the  Moore- Penrose  pseu¬ 
doinverse.  Finally,  the  pole  estimates  are  found  by 

Pj  m  •  i  =  1,1,. -..A-  (7) 

where  P(r)  s:  1  +  4-  •••  +  lit1.  Once  the  poles  have 

been  determined,  tbe  amplitude  coefficients  can  be  found 
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nbiaaad  aatimaton  rise*  it  ia  Uaad  as  tbc  CRB.  The 
resolution  limit  lor  the  TLS-Prae;  aatimatioo  algorithm 
it  ■imiUrly  pro:  by  ariag  pole  aagk  variasee  statistics 
lor  the  tne  poles;  these  are  readily  obtained  from  o  -  Z§ 
(see  (»)). 


Bound  An 


sis  for  Tore  Vadtap 


Note  that  £  mode  estimates  are  obtained,  o i  which  n 
are  "true*  modes.  For  £  >  n,  £  -  n  of  the  mode  esti¬ 
mates  are  extraneous.  If  n  is  knows  then  the  true  modes 
can  be  identified  at  the  n  highest  energy  estimates  (3). 
Howeeer,  ia  practice  n  is  typically  unknown;  ia  this  cate 
the  number  of  singular  values  retained  it  a  fixed  upper 
bound  of  n,  and  true  modes  could  be  separated  from  ex¬ 
traneous  modes  by  using  a  mode  energy  threshold  (as 
discussed  below). 


In  order  to  establish  a  resolution  bound  and  probabili¬ 
ties  of  detection  and  false  alarm,  we  need  the  statistics 
for  the  estimated  parameters  in  the  TLS-Prony  model. 
For  the  resolution  bounds  and  detection  probabilities,  we 
need  only  the  statistics  for  the  "true’'  modes  (that  it,  the 
n  modes  srith  highest  energies).  For  this  case,  the  statis¬ 
tics  of  the  parameter  estimates  have  been  found  in  (3,7). 
However,  for  the  false  altsrm  probability,  we  also  need 
the  statistics  of  the  £  — a  extraneous  modes.  The  deriva¬ 
tion  of  the  extraneous  mode  statistics  appears  in  (8);  the 
main  result  is  summarised  in  the  following  Theorem. 

Theorem  1:  Assume  y  is  as  given  in  Equation  1.  Let  p  = 
{p.  and  i  sr  be  at  in  Equation  1,  and  let  p'  = 

{p* )  and  x*  m  {*  J  }^~*  denote  the  £  -  a  extraneous 
modes  obtained  is  the  TLS-Prony  procedure  when  0  s 
0.  Let  0  denote  the  4£  X  1  vector  containing  the  angles 
and  ^magnitudes  of  p,  x,  p\  and  respectively,  and 
let  0  denote  the  TLS-Prony  estimate  of  I.  Then  the 
asymptotic  (as  o  — •  0)  probability  density  function  of  0 
is  Gaussian:  ^ 

d~N(«,o-E,),  (9) 

where  E#  is  a  covariance  matrix  which  depends  on  m,  £, 

•»<! 

The  proof  of  tbc  Theorem  ud  as  explicit  expression  for 
£#,  can  be  found  in  (8).  Also  in  (8)  is  an  expression  for 
£#  when  €  is  reparameteriaed  in  term*  of  the  real  and 
imaginary  part*  of  the  poles  and  amplitude  coefficient*. 
We  note  that  the  above  Theorem  give*  a  theoretical  ex¬ 
pression  for  the  complete  pdf  of  the  estimated  parame¬ 
ter*;  this  pdf  can  then  be  used  to  study  the  resolution 
properties  and  detection  probabilities  of  scattering  cen¬ 
ters,  as  is  discussed  below. 


ANAIYSIS 


Consider  two  poles  on  the  unit  circle,  p,  b  oj  r'"'1  and 
Pi  *  oj  We  define  the  resolution  limit  r  as 

r  se  2ow,  <+2*7*,,  (10) 

where  the  angle  variances  of  pi  and  pa  are  *1,  and  **  , 
respectively.  When  two  poles  are  at  this  limit,  the  9b% 
confidence  intervals  of  the  angle  estimates  for  each  pole 
become  disjoint. 

The  CRB  resolution  limit  is  found  by  using  CRB  ex¬ 
pressions  for  the  pole  angles  in  two  pole  model.  The 
CRBs  for  the  pole  angles  are  inserted  into  oWi  and 
in  Equation  10;  such  expressions  can  be  found,  for  ex¬ 
ample,  in  (0).  This  limit  gives  a  lower  bound  for  all 


la  this  study  we  consider  two  aqua)  energy  modes  located 
os  the  unit  circle  at  pj  *  M(]  p,  B 

lor  a  data  length  of  m  *  10,  where  /  Is  the  separation 
of  the  two  modes  in  Fbcrier  bins.  Fbr  the  TLS-Prony 
simulations  a  prediction  order  of  L  *  4  was  need. 

Figure  1  shows  the  bounds  far  the  CRB,  TLS-Prony  sta¬ 
tistical  theory,  aad  TLS-Prony  Monte-Carlo  results.  The 
axes  are  normalised  to  make  the  corves  independent  of 
data  length  fbr  the  CRB  curve  aad  the  TLS-Prony  sta¬ 
tistical  theory  curve.  From  these  corves  we  can  see  that 
the  TLS-Prony  theoretical  bound  is  quite  close  to  the 
CRB  over  a  wide  range  of  SNR.  Recall  that  the  CRB 
and  TLS-Prony  bounds  are  both  derived  using  small 
perturbation  analysis,  so  hold  only  fbr  high  SNR.  The 
TLS-Prony  Monte  Carlo  simulations  show  good  agree¬ 
ment  with  the  theoretical  bound  above  18  dB  SNR.  We 
can  see  that  above  20dB  SNR/pole/bin  the  TLS-Prony 
method  virtually  achieves  the  CRB. 

Below  20dB  SNR/pole /bin,  the  Monte-Carlo  simulations 

?‘ve  much  higher  variance  than  both  the  theoretical  TLS 
rosy  curve  and  the  CRB  curve.  Note,  however,  that 
the  TLS-Prony  analytical  variance  expression  was  de¬ 
rived  under  the  assumption  of  high  SNR,  aad  is  not  ex¬ 
pected  to  be  accurate  at  low  SNR.  In  addition,  the  CRB 
is  not  necessarily  a  tight  bound  at  km  SNR.  Thus,  it 
is  not  dear  what  the  minimum  achievable  variance  is  in 
this  region.  For  example,  it  is  not  known  whether  (or 
how  much)  an  iterative  maximum  likelihood  procedure 
would  result  ia  lower  variance  in  this  region. 

We  note  that  above  about  18  dB  SNR/pole/bin,  the 
TLS-Prony  technique  gives  resolutions  of  less  than  one 
Fourier  bin.  Therefore,  the  resolution  of  the  TLS-Prony 
technique  is  better  than  FFT-based  techniques  since 
P FT- based  techniques  can  only  resolve  to  within  one 
Fourier  bin.  If  windowing  is  used  in  conjunction  with 
the  FFT-based  methods,  their  resolution  is  even  larger 
than  one  Fourier  bin  (e.f.,  it  would  be  about  1 J  Fourier 
bins  using  a  Bamming  window). 


In  practice,  ose  does  not  know  s  priori  how  many  scat¬ 
tering  centers  are  present.  In  this  case,  one  would  accept 
or  reject  a  mode  estimate  as  a  scattering  center  based  on 
some  threshold.  We  consider  a  threshold  on  the  energy 
of  the  estimated  exponential  mode,  as  this  corresponds 
to  radar  crass  section  of  an  estimated  scattering  center. 

In  the  TLS-Prony  method,  one  obtains  estimates  of  L 
poles  and  L  corresponding  amplitude  coefficients.  FVom 
this,  one  can  compute  the  energy  £t  of  each  of  the  L 
modes  by 

;■=  j.j . (ii) 


where  0,  and  a,  are  the  magnitudes  of  the  jth  ampli¬ 
tude  coefficient  and  pole,  respectively.  We  consider  an 
estimated  mode  to  be  detected  as  a  valid  scattering  cen¬ 
ter  if  it*  energy  exceeds  a  prespecified  threshold,  aad  we 
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reject  the  mode  m  an  invalid  scattering  center  if  it  doe* 
not.  We  then  can  present  detection  result*  is  the  form 
a i  receiver  operation  characteristic  (ROC)  curves. 

Probability  of  Detection 

We  define  a  detection  to  he  the  caae  in  which  all  of  the 
true  mode  energy  estimates  exceed  an  energy  threshold, 
E°.  We  thus  now  derive  the  energy  statistic*  for  the 
true  modes  (i.e.,  j  =  1,2,...,*).  These  statistics  can 
be  found  from  the  statistical  pdf  given  in  Theorem  1. 
It  can  be  shown  (6)  that  the  energies  are  noo central  x* 
distributed;  for  high  SNR,  the  aoacentral  x*  distribution 
is  well-approximated  by  a  Gaussian  distribution.  Using 
first-order  approximations,  it  is  possible  to  derive  the 
mean  and  covariance  of  this  distribution.  Thus,  from 
Theorem  1  we  have  the  following  corollary. 

Corollary  J:  Let 

S=|£,  E,  ■■■  E.f,  (12) 

denote  the  parameter  vector  for  the  mode  energies  of 
the  true  modes  (s.e.,  the  mode  energies  corresponding 
to  p  and  x).  Let  E  denote  the  estimated  energies  corre¬ 
sponding  to  the  TLS-Prony  parameter  estimates.  Then 
the  asymptotic  (as  o  — »  0)  pdf  of  £  is  given  by 

2~N(E,e  Er),  (IS) 

where  £c  depends  on  m,  L ,  and  (x,,p,}*Ml.  An  explicit 
expression  for  E;  can  be  found  in  (8). 

Given  the  true  mode  energy  distribution,  the  probabilty 
of  detection,  Pd,  it  given  by 

rc  =  p  •j(ti>Ec,E,>E° . £>s“).  (u) 

This  probability  is  readily  computed  using  Equation  13. 

To  verify  that  the  theoretical  energy  distribution  given 
above,  Monte- Carlo  simulations  were  performed  for  a 
two  mode  case.  In  this  case,  the  data  consists  of  two 
equal  energy  modes,  with  ij  =  *j  =  1,  located  on  the 
unit  circle  spaced  one  Fourier  bin  apart  at  px  s= 
and  pj  =  foy  i  data  length  of  m  =  10.  A  pre¬ 

diction  order  of  I  c  4  was  used,  and  two  singular  val¬ 
ues  were  retained.  The  SNR  for  these  simulations  was 
lOdB/pole. 

Figure  2  shows  a  comparison  between  the  theoretical  pdf 
and  a  histogram  obtained  from  Monte-Carlo  simulations 
(note  that  both  mode*  have  the  same  theoretical  pdfs 
and  had  similar  histograms).  It  can  be  seen  that  the 
theoretical  energy  distribution  is  a  good  approximation 
to  simulation  results  in  this  case. 

Probability  ef  Fblse  Alarm 

We  define  a  false  alarm  to  be  the  case  in  which  one  or 
more  of  the  extraneous  modes  is  above  the  energy  thresh¬ 
old,  E°%  and  thus  misidentified  as  a  true  mode.  We 
can  derive  the  statistal  properties  of  the  estimated  ex¬ 
traneous  mode  energies  in  a  similar  manner  as  above.  In 
this  case,  however,  the  “true*  energies  of  the  extraneous 
modes  are  tero,  so  the  extraneous  modes  are  distributed 
as  centra]  Chi-squared  with  two  degree*  of  freedom,  xl 
(see  (8)).  This  is  stated  in  the  following  corollary. 

Corollary  t:  Let  j  denote  the  estimated  en¬ 
ergies  corresponding  to  the  TLS-Prony  parameter  esti¬ 


mates  of  the  extraneous  mode  energies.  Then  the  asymp¬ 
totic  (a*  o  -•  0)  pdfs  of  these  energies  are  given  by 

£*  ~  *5  (* '  E*;)  s  *  1,2,...,!-  n.  (15) 

where  £r  does  not  depend  on  o.  An  explicit  expression 
lor  £*•  can  be  found  in  (6). 

Given  the  extraneous  mode  energy  distributions,  the 
probabilty  of  false  alarm,  Pta%  b  given  by 

Eta  « 1  -  Pi  (£}  <  £°,£J  <  ^ . C.  <  **)  • 

06) 

Note  that  since  the  extraneous  mode  energy  distribution 
is  Chi-squared  with  two  degrees  of  freedom,  Pr  a  can  be 
evaluated  using  a  Rayleigh  distribution. 

Figure  3  shows  a  comparison  betwom  the  theoretical  pdf 
lor  the  extraneous  modes  in  the  previous  two  modes 
example  and  a  histogram  obtained  from  Monte-Carlo 
simulations.  Note  that  the  theoretical  predictions  agree 
closely  with  Monte-Carlo  simulations. 

ROC  Analysis  for  Two  Undamped  Modes 

Using  the  above  detection  and  false  alarm  probability 
results,  we  can  derive  ROC  curves  lor  scattering  cen¬ 
ter  detection  at  various  SNR*.  Figure  4  presents  such 
curves  for  the  case  considered  in  Figures  2  and  3,  but 
with  varying  SNR  (note  that  1  —  Pp  is  actually  plotted 
along  the  vertical  axis).  For  an  SNR  at  or  above  lOdB 
per  pole,  Pd  is  always  above  0.9  even  when  Pta  is  very 
small  ( t.g.y  1C7).  However  this  is  not  the  case  for  lower 
SNR.  Note  that  for  low  SNRs  Pd  never  reaches  one  even 
if  Pta  is  one.  This  is  because  the  true  mode  energy  dis¬ 
tributions  were  approximated  as  Gaussian,  and  the  tail 
of  this  Gaussian  distribution  gives  a  nonzero  probability 
of  a  negative  energy  (the  noncentral  x\  distribution  does 
not  have  such  a  tail).  For  high  SNR,  the  approximation 
becomes  more  valid. 

In  computing  the  curves  in  Figure  4,  it  was  assumed 
that  the  extraneous  mode  energy  distributions  are  inde¬ 
pendent.  Note  that  this  is  a  worst  case  assumption,  since 
Pta  would  decrease  if  the  extraneous  mode  energies  were 
dependent. 

CONCLUSIONS 

In  the  paper  we  presented  resolution  hounds  and  detec¬ 
tion  results  for  estimating  the  scattering  centers.  These 
bounds  are  based  on  an  exponential  model  of  target  scat¬ 
tering.  which  generalizes  the  point  scattering  model.  The 
popular  TLS-Prony  algorithm  was  used  to  estimate  the 
parameters  of  the  exponential  model.  A  high  SNR  sta¬ 
tistical  analysis  of  the  TLS-Prony  algorithm  was  first 
presented.  Then,  based  on  the  results  of  the  statistical 
analysis,  both  resolution  bounds  and  detection  results 
were  presented.  The  resolution  bounds  were  compared 
with  both  the  Cramer- Rao  Bound  and  with  Monte-Carlo 
simulations.  It  was  shown  that  for  an  SNR  above  18  dB, 
the  TLS-Prony  method  is  capable  of  resolution  to  less 
than  a  Fourier  bin. 

The  probabilities  of  the  detection  and  the  the  false  alarm 
were  derived  based  on  the  mode  energy  distributions  for 
both  the  true  and  the  extraneous  modes.  For  high  SNR 
true  mode  energy  distributions  can  be  approximated  as 
Gaussian  distributions,  but  extraneous  mode  energy  dis¬ 
tributions  are  central  chi-squared  distributed.  These  de¬ 
tection  and  false  alarm  probabilities  can  be  presented  as 
ROC  curvet  for  scattering  center  detection  for  examples 
of  interest. 
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Figure  1:  Resolution  bounds  for  two  equal  energy  un¬ 
damped  modes. 


Figure  2:  Triie  mode  theoretical  pdf  and  histogram  for 
two  equal  energy  undamped  modes. 


Figure  3:  Extraneous  mode  theoretical  pdf  and  his¬ 
togram  for  two  equal  energy  undamped  modes. 


Figure  4:  ROC  curves  for  two  equal  energy  undamped 
modes. 
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D.  Reprint  of  “Prony  Modeling  of  Linear  FM  Radar  Data” 


The  following  pages  contain  a  preprint  of  the  technical  report,  “Prony  Modeling  of 
Linear  FM  Radar  Data”. 
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Abstract 


In  this  report  some  preliminary  research  results  on  the  parametric  modeling  of  linear 
FM  radar  data  are  presented.  Superimposed  exponential  signals  are  used  to  model 
the  data.  The  model  parameters  are  then  estimated  by  the  TLS-Prony  estimator.  In 
addition  to  the  TLS-Prony  algorithm,  a  modified  algorithm  incorporating  the  TLS- 
Prony  estimator  and  a  residue  concept  is  proposed  to  estimate  the  parameters.  The 
estimates  from  the  modified  algorithm  are  then  improved  by  a  maximum  likelihood 
estimator  since  the  residue  algorithm  introduces  bias  on  the  estimates.  The  estimation 
results  are  compared  to  the  results  from  the  fast  Fourier  transform  technique.  The 
issue  of  reduced  bandwidth  is  also  considered  in  this  report. 


105 


I.  Introduction 


This  report  presents  some  preliminary  research  on  the  parametric  modeling  of 
linear  FM  radar  data.  Based  on  the  scattering  center  assumption  the  data  are  modeled 
by  superimposed  exponential  signals  [1], 

n 

yq -^2xiPql  +  eg  q  =  0,1,. ..,m  —  1,  (1) 

i= 1 

where  n  is  the  number  of  the  scattering  centers  and  m  is  the  number  of  frequency 
measurements.  In  addition,  each  p,  is  a  pole  corresponding  to  a  scattering  center, 
each  a,  is  the  amplitude  associated  with  that  pole,  and  eq  is  measurement  noise. 

The  problem  is  to  estimate  the  parameters  in  the  exponential  model.  Many  al¬ 
gorithms  have  been  developed  to  achieve  this  objective.  One  of  the  most  popular 
methods  is  the  TLS-Prony  algorithm  [2].  Numerical  studies  and  statistical  analysis 
have  shown  that  for  high  SNR,  the  TLS-Prony  estimator  is  unbiased  and  has  very 
promising  performance.  Details  about  the  TLS-Prony  analysis  can  be  found  in  [3]. 
In  this  report  we  are  concerned  about  how  the  TLS-Prony  algorithm  performs  with 
respect  to  the  traditional  fast  Fourier  transform  (FFT)  technique,  when  applied  to 
radar  data.  For  this  study  we  use  linear  FM  radar  data  of  an  aircraft  that  was 
collected  at  Rome  Laboratories. 

In  addition,  a  modified  algorithm  incorporating  TLS-Prony  and  a  residue  concept 
is  proposed.  The  parameter  estimates  from  the  modified  algorithm  then  are  improved 
by  a  maximum  likelihood  estimator  (MLE)  since  the  residue  algorithm  introduces  bias 
on  the  parameter  estimates.  The  MLE  presented  here  is  a  complex  version  of  the  MLE 
presented  in  [4]. 
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Bandwidth  reduction  is  also  considered  in  the  report.  Theoretically,  parametric 
modeling  has  no  resolution  limit  (practically,  of  course,  this  is  not  the  case).  Here 
we  would  like  to  know  how  parametric  modeling  performance  degrades  as  bandwidth 
is  reduced.  The  TLS-Prony  estimation  performance  with  reduced  bandwidth  is  com¬ 
pared  to  full  bandwidth  performance  of  the  FFT  technique. 

An  outline  of  this  report  is  as  follows.  Section  II  presents  the  TLS-Prony  algorithm 
and  the  modified  algorithm.  In  Section  III,  the  estimated  results  from  both  algorithms 
are  compared  to  each  other  and  to  the  FFT  technique  as  well.  Section  IV  concludes 
this  report. 


II.  Estimation  Algorithms 


A.  TLS-Prony  Algorithm 

Based  on  the  data  sequence  in  Equation  (1),  backward  linear  prediction  equations 
are  formulated  as: 


where 


(2) 


(3) 
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and  where 


t/o  2/i  2/2  •  •  •  2 IL 

2/i  2/2  2/3  •  •  •  2/L+i 

2/m-(L+l)  Vm—L  Vm-iL-l)  '  '  ’  2/m-l 

Here  L  is  the  order  of  prediction  and  b  is  the  coefficient  vector  of  the  polynomial  B(z) 
given  by 

B{z)  =  l  +  blz  +  b2z2  +  ---  +  bLzL.  (5) 

For  the  noiseless  case,  L  can  be  any  integer  greater  than  or  equal  to  the  model  order 
n;  however,  choosing  L  >  n  results  in  more  accurate  parameter  estimates. 

The  solution  of  Equation  (2)  involves  obtaining  an  SVD  of  the  matrix  y  :  Y 
and  truncating  all  but  the  first  n  singular  values  to  arrive  at  an  estimate  y  •  Y 
[2],  This  leads  to  the  modified  linear  prediction  equation 

Yb  =  -y  (6) 

from  which  the  linear  prediction  coefficient  vector  estimate  b  is  found  as 

b=-Y+y ,  (7) 

where  +  denotes  the  Moore-Penrose  pseudoinverse.  Finally,  the  estimates  for  the 
poles  are  found  by 

Pj  =  zeroj  ,  j  =  1,2,...,  L.  (8) 

Note  that  for  high  frequency  radar  target  responses,  the  scattering  is  specular  and 
thus  the  poles  are  very  near  the  unit  circle.  That  is,  the  magnitudes  of  the  estimated 
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poles  are  close  to  unity.  Before  estimating  the  amplitude  coefficients  we  thus  force 
the  estimated  poles  to  be  on  the  unit  circle  (by  setting  the  magnitude  to  be  one).  It 
is  believed  this  constraint  on  the  poles  may  give  better  results.  Here  we  will  focus  on 
the  constraint  case. 

Once  the  poles  have  been  determined,  the  amplitude  coefficients  can  be  found 
using  the  pole  estimates  and  Equation  (1).  This  leads  to  the  following  least  squares 
equation  for  the  amplitude  coefficients, 


or 

AlXl  ~  Ya.  (10) 

The  amplitude  coefficients  can  be  found  from  a  least  squares  solution  to  Equation  (10), 

XL  =  (aial)-' A-Ly.  =  AtY.,  (II) 

where  ’  denotes  complex  conjugate  transpose.  Because  only  n  singular  values  of 
y  •  Y  are  nonzero,  there  are  at  most  n  pole  estimates  which  can  correspond  to 
data  modes.  Therefore,  only  the  n  poles  which  have  the  largest  energy  are  retained. 
This  is  done  by  computing  the  L  mode  energies  as 

N  m  —  1 

£.  =£l*.f  £!p.|2'  *=  1,2 . i  (12) 

t= 1  q= 0 
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and  retaining  those  n  poles  whose  corresponding  energies  are  highest.  We  then  re- 
estimate  the  amplitude  coefficients  of  these  n  poles.  This  is  done  using  Equation  (11), 
except  that  AL  is  replaced  by  A,  where  A  is  the  Vandermonde  matrix  composed  only 
of  the  n  columns  of  Al  corresponding  to  the  n  highest  energy  poles. 

B.  A  Modified  TLS-Prony  Algorithm  Incorporating  a  Residue  Concept 

A  modified  algorithm  incorporating  the  TLS-Prony  algorithm  with  a  residue  con¬ 
cept  is  proposed  here.  After  estimating  the  parameters  via  the  TLS-Prony  algorithm, 
we  can  use  Equation  (1)  to  create  the  estimated  data  sequence.  Assuming  perfect 
modeling  on  the  true  signals  the  difference  between  the  measurement  data  sequence 
and  the  estimated  data  sequence  would  be  the  noise  sequence  only.  However,  since 
the  modeling  is  not  perfect,  the  residue  contains  some  information  about  the  data. 
This  motivates  the  following  algorithm. 

•  Lise  the  estimated  parameters  from  TLS-Prony  to  create  the  estimated  data 
sequence  via  Equation  (1).  Subtract  this  sequence  from  the  original  sequence, 
giving  a  residual  data  sequence. 

•  Use  the  TLS-Prony  algorithm  on  the  residual  data  sequence  to  obtain  residue 
model  parameters.  Combine  the  two  sets  of  parameters  (the  set  fiom  the  original 
measurement  data,  and  the  set  from  the  residual  data). 

•  Use  a  maximum  likelihood  estimator  (MLE)  to  improve  the  combined  estimates. 

The  MLE  is  used  to  overcome  bias  resulting  from  estimating  the  two  parameter 
sets  separately.  The  MLE  chosen  is  the  complex  version  of  the  MLE  presented  in  [4]. 
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III.  Results 


In  this  report  a  Hamming  window  is  always  used  to  reduced  the  sidelobes  at 
the  cost  of  losing  resolution.  The  results  shown  here  will  only  focus  on  one  of  the 
seventeen  data  sequences  since  these  sequences  are  responses  for  the  same  target,  and 
are  nearly  identical  to  one  another. 

A.  Some  Comparisons  using  the  data  sequence  wbOOOl 

In  Figure  1  the  magnitude  of  the  inverse  FFT  (IFFT)  of  the  1024  point  data 
sequence  wbOOOl  using  a  Hamming  window  is  shown.  The  high  energy  section  of  this 
IFFT  is  shown  in  Figure  2.  We  will  use  this  result  as  a  standard  to  compare  the 
estimates  from  both  the  constrained  TLS-Prony  algorithm  and  the  modified  TLS- 
Prony  coupled  with  the  residue  and  MLE  algorithm. 

The  estimation  result  for  the  constrained  TLS-Prony  algorithm  using  the  original 
data  is  compared  to  the  IFFT  of  the  original  data  in  Figure  3  (the  solid  line  is  for  the 
estimated  data  and  the  broken  line  is  the  original  data).  Although  the  estimated  data 
do  not  fit  the  original  data  very  well,  from  Figure  4,  which  shows  the  line  spectrum  of 
the  TLS-Prony  estimates,  we  can  see  that  at  some  points  (as  indicated  in  the  figure) 
we  superresolve  closely  spaced  scattering  centers. 

For  the  modified  algorithm  the  results  are  even  better  and  shown  in  Figures  5 
and  6.  In  Figure  5  it  is  clear  that  the  estimated  data  (solid  line)  fit  the  original  data 
(broken  line)  very  well.  The  line  spectrum  in  Figure  6  shows  the  superresolvability 
of  the  modified  algorithm.  Many  scattering  centers  not  seen  in  the  original  data  are 
clear  now.  This  extra  resolution  is  obtained  at  the  cost  of  increased  computation. 
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In  the  modified  algorithm  we  use  TLS-Prony  twice,  which  doubles  the  computation. 
In  addition,  the  computational  cost  of  the  MLE  has  to  be  added.  On  average  the 
computation  for  the  modified  algorithm  is  four  times  that  of  the  computation  for  the 
constrained  TLS-Prony  algorithm. 
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Figure  2 
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Figure  1:  The  IFFT  of  the  original  data  sequence. 


The  IFFT  of  the  original  data  sequence  wbOOOl  (zoom) 


:  The  high  energy  section  of  the  IFFT  of  the  original  data  sequence. 
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Figure  3:  The  IFFTs  of  the  original  data  sequence  and  the  estimated  data  sequence 
by  the  TLS-Prony  algorithm  using  the  original  data. 


Figure  4:  The  line  spectrum  of  the  TLS-Prony  estimates  using  the  original  data. 
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Figure  5:  The  IFFTs  of  the  original  data  sequence  and  the  estimated  data  sequence 
by  the  modified  algorithm  using  the  original  data. 

The  line  spectrum  of  the  estimated  sequence  (mod.  algo.  (1024)) 
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Figure  6:  The  line  spectrum  of  the  modified  algorithm  estimates  using  the  original 
data. 
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B.  Reduced  bandwidth  consideration 


In  the  previous  subsection  we  have  shown  that  the  modeling  techniques  based  on 
the  TLS-Prony  algorithm  give  very  promising  results,  and  outperform  the  conven¬ 
tional  FFT  technique  in  terms  of  resolving  closely  spaced  scattering  centers.  This 
motivates  the  following  question.  If  a  smaller  set  of  the  original  data  sequence  is 
used,  is  the  performance  of  these  modeling  techniques  still  good  enough  to  compete 
with  the  FFT  techniques  using  full  bandwidth?  We  investigate  this  question  using 
some  results  based  on  the  same  set  of  data  as  before. 

In  Figure  7  the  IFFTs  of  the  original  data  (broken  line)  and  the  estimates  from 
the  TLS-Prony  algorithm  using  the  middle  512  points  (solid  line)  are  shown.  The 
corresponding  line  spectrum  is  shown  below  in  Figure  8.  Here  we  can  see  that  the 
estimate  is  not  impressive.  In  fact,  the  TLS-Prony  algorithm  misses  some  of  the 
scattering  centers  (as  indicated  in  the  figure).  However  the  results  from  the  modified 
algorithm  are  much  better.  They  are  shown  in  Figures  9  and  10.  These  figures  are 
analogous  to  Figures  7  and  8  except  the  modified  algorithm  is  used.  It  is  clear  from 
these  figures  that  the  modified  algorithm  obtains  all  the  scattering  centers  shown 
for  the  FFT  technique,  and  it  also  even  superresolves  some  closely  spaced  scattering 
centers.  This  suggests  that  the  required  bandwidth  can  be  reduced  by  a  factor  of  two 
when  the  modified  algorithm  is  used. 
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IFFTs  of  the  estimated  sequence  (TLS-Prony(512))  and  the  data  (1024) 
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Figure  7:  The  IFFTs  of  the  original  data  sequence  and  the  estimated  data  sequence 
by  the  TLS-Prony  algorithm  using  the  middle  512  points  of  the  original  data  . 
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Figure  8:  The  line  spectrum  of  the  TLS-Prony  estimates  using  the  middle  512  points 
of  the  original  data. 
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EFFTs  of  the  estimated  sequence  (mod.  alg.(512))  and  the  data  (1024) 
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Figure  9:  The  IFFTs  of  the  original  data  sequence  and  the  estimated  data  sequence 
by  the  modified  algorithm  using  the  middle  512  points  of  the  original  data. 


The  line  spectrum  of  the  estimated  sequence  (mod.  algo.  (512)) 
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Figure  10:  The  line  spectrum  of  the  modified  algorithm  estimates  using  the  middle 
512  points  of  the  original  data. 
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IFFTs  of  the  estimated  sequence  (TLS-Prony(256))  and  the  data  (1024) 


Figure  11:  The  IFFTs  of  the  original  data  sequence  and  the  estimated  data  sequence 
by  the  TLS-Prony  algorithm  using  the  middle  256  points  of  the  original  data. 

In  Figures  11,  12,  13,  and  14  we  show  the  results  for  using  the  middle  256  points 
of  the  original  data  sequence.  For  the  TLS-Prony  algorithm  (in  Figures  11  and  12) 
several  scattering  centers  are  missed.  For  the  modified  algorithm,  however,  the  results 
are  more  encouraging.  These  results  are  shown  in  Figures  13  and  14.  Some  scattering 
centers  are  still  missed,  but  not  the  stronger  scattering  centers. 
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The  line  spectrum  of  the  estimated  sequence  (TLS-Prony  (256)) 


Normalized  Range 


Figure  12:  The  line  spectrum  of  the  TLS-Prony  estimates  using  the  middle  256  points 
of  the  original  data. 


lFPTs  of  the  estimated  sequence  (mod. algo. (256))  and  the  data  (1024) 
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Figure  13:  The  IFFTs  of  the  original  data  sequence  and  the  estimated  data  sequence 
by  the  modified  algorithm  using  the  256  points  of  the  original  data. 
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IV.  Conclusions 


In  this  report  we  presented  some  preliminary  results  on  parametric  modeling  of 
linear  FM  radar  data  collected  at  Rome  Laboratories.  Modeling  techniques  and  the 
conventional  FFT  technique  were  compared.  We  focused  on  the  TLS-Prony  algorithm 
and  proposed  a  modified  algorithm  incorporating  TLS-Prony,  a  residue  concept,  and 
MLE.  Several  issues  were  discussed  in  the  report.  First,  for  the  same  bandwidth, 
the  modeling  techniques  used  in  general  outperform  the  FFT  technique,  especially 
the  modified  algorithm.  In  addition,  we  found  that  using  the  modified  algorithm 
we  may  reduce  the  bandwidth  by  a  factor  of  two  without  losing  the  performance  in 
comparison  to  the  FFT  technique.  Further  investigation  is  needed  on  the  modified 
algorithm  in  order  to  fully  understand  its  performance. 
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